• Self-consistent mode-coupling theory for the viscosity of ...

Self-consistent mode-coupling theory for the viscosity of rodlike
polyelectrolyte solutions
Kunimasa Miyazakia)
Department of Chemistry, Columbia University, New York, New York 10027
Biman Bagchi
Solid State and Structural Chemistry Unit, Indian Institute of Science, Bangalore 560 012, India
Arun Yethirajb)
Theoretical Chemistry Institute and Department of Chemistry, University of Wisconsin, Madison,
Wisconsin 53706
A self-consistent mode-coupling theory is presented for the viscosity of solutions of charged rodlike
polymers. The static structure factor used in the theory is obtained from polymer integral equation
¨
theory; the Debye-Huckel approximation is inadequate even at low concentrations. The theory
predicts a nonmonotonic dependence of the reduced excess viscosity R on concentration from the
behavior of the static structure factor in polyelectrolyte solutions. The theory predicts that the peak
in R occurs at concentrations slightly lower than the overlap threshold concentration, c * . The peak
height increases dramatically with increasing molecular weight and decreases with increased
concentrations of added salt. The position of the peak, as a function of concentration divided by
c * , is independent of salt concentration or molecular weight. The predictions can be tested
experimentally.
I. INTRODUCTION and Schmidt,1 this problem has witnessed some grave mis-
takes, causing much confusion.
Polyelectrolyte solutions are widely considered to be one A contributing factor to the above-mentioned confusion
of the least understood substances in polymer science.1 There is the fact that a review of experimental data on the viscosity
are several features of these solutions that make them rather of polyelectrolyte solutions shows large inconsistencies.4,5 It
complex. For one, the long-ranged nature of the electrostatic is only recently that the reason for this discrepancy between
interactions results in long-ranged correlations even in dilute different experimental measurements has been established.
solutions. In addition, polymer conformations are very sen- In careful measurements of the shear rate dependence of the
sitive to concentration and ionic strength because the elec- viscosity, Boris and Colby4 and Krause et al.5 showed that
trostatic interactions compete with short-ranged ‘‘hydropho- polyelectrolyte solutions were shear thinning at extremely
bic’’ interactions. This complex static behavior is low shear rates, and argued that most of the older experi-
accompanied by very interesting dynamical behavior of ments did not report the relevant viscosity in the low shear
polyelectrolyte solutions. In this work, we present a theoret- limit.
ical study of the viscosity of dilute and semidilute rodlike In the last ten years, several theories studies have been
polyelectrolyte solutions using a liquid state approach. put forward. Notable among them are the mode coupling
The viscosity of polyelectrolyte solutions displays an theory calculation by Borsali et al.6 and by Rabin et al.,3 the
interesting ‘‘anomalous’’ concentration dependence at scaling theory of Dobrynin et al.,7 the effective medium
low polymer concentrations.2,3 The quantity that is normally theory of Muthukumar,8 and the Kirkwood theory of Nishida
discussed is the reduced viscosity R defined as R et al.9 The scaling theory does not predict a peak in the re-
( 0 )/( 0 c p ), where is the viscosity of the solution, duced viscosity, in salt-free solutions. The theory of Muthu-
0 is the viscosity of the solvent in the absence of the poly- kumar argues that the peak in R arises from a screening of
mer, and c p is the monomer concentration of polymers. intramolecular hydrodynamic interactions as the concentra-
Experiments show that R displays a sharp peak at low poly- tion is increased. The theories of Borsali et al.,6 Rabin
electrolyte concentrations for a variety of different solutions. et al.,3 and Nishida et al.9 predict that the peak in R arises
This anomalous concentration dependence of viscosity of di- due to increased screening from counterions as the concen-
lute polyelectrolyte solutions has been the focus of attention tration is increased, and are similar in spirit to the mode-
for over 50 years and, although there have been many theo- coupling theory for charged colloids10 which is argued to be
ries that address the problem, it is not considered to be well accurate for the viscosity of spherical polyelectrolytes.11
understood.1,4,5 In fact, as discussed eloquently by Forster Other theories rely on conformational changes of polymers
with changing concentration. These theories do not take into
a
Electronic mail: km2233@columbia.edu account the interesting behavior in the static structure factor
b
Electronic mail: yethiraj@chem.wisc.edu in dilute polyelectrolyte solutions.
 The basic idea of the present work is that the features in The rest of this paper is organized as follows: Section II
R arise from the behavior of the static structure factor of outlines the theory, Sec. III presents some results and a dis-
dilute and semidilute polyelectrolyte solutions. In dilute so- cussion, and Sec. IV presents some conclusions.
lutions, the static structure factor displays a prominent peak
at low wave vectors.12,13 As the concentration is increased II. MODE-COUPLING THEORY
the peak broadens and moves to higher wave vectors see
Fig. 2 and discussion in Sec. III . This indicates the presence The polyelectrolyte solution consists of charged rodlike
of strong liquid like correlations on long length scales at low polymers and counterions. Each rod consists of N tangent
concentrations; correlations that become less important as charged hard spheres of diameter ; counterions are charged
the concentration is increased. This observation naturally hard spheres of diameter . We combine a recently devel-
leads to the question: Could this strong non monotonic con- oped quantitatively accurate theory for the liquid structure in
centration dependence of static pair correlations be the main charged interacting polyelectrolytes18 –20 with a self-
physics behind the anomalous behavior of the viscosity? To consistent mode coupling theory MCT to study the dynam-
this end, we develop a liquid state theory that incorporates ics. The starting point for the calculation of viscosity is the
the behavior of the static structure factor of polyelectrolyte Green-Kubo formula21
solutions. While it has been suggested14,15 that the interesting 1 zx zx
concentration dependence of viscosity in dilute polyelectro- lim dt k,t k,0 , 2.1
k→0 k BTV 0
lyte solutions could arise from the intermolecular pair corre-
lations as reflected in the peak in the structure factor at where k B is Boltzmann’s constant, T is the temperature, V is
small wave numbers , the relationship between the viscosity the volume, zx (k,t) is the transverse or off-diagonal com-
and the structure factor is by no means obvious. It is of ponent of the wave vector k and time t dependent stress
interest, therefore, to develop a quantitative theory relating tensor, and ¯) denotes an average over an equilibrium en-
the static structure to the viscosity. semble. The total transverse stress tensor of a polyelectrolyte
We consider a system of charged rods and present a self- solution contains contributions from solvent, polymer, and
consistent mode-coupling theory for the viscosity. We choose small ions.
to study rods in order to focus on intermolecular effects. For dilute and semidilute polyelectrolyte solutions sev-
There have been many theories that explain the concentra- eral simplifications are possible: First of all, the contribution
tion dependence of the reduced viscosity on intramolecular of the solvent is simply given by the viscosity in the absence
effects. Studying a system with rigid molecules allows us to of the solute 0 . It is because the presence of low concen-
isolate intermolecular effects since conformational changes trations of polymer and electrolyte are not expected to alter
and intramolecular interactions are absent. Although we are the solvent dynamics. Second, there is a contribution of the
not aware of experiments for the viscosity of rodlike poly- rotational Brownian motion.22 Last, there is a contribution
electrolytes, these are certainly possible, for example on so- due to polymer-polymer interactions p-p which is expected
lutions of tobacco mosaic virus TMV particles. to dominate over the contributions from small ions. There-
We present a self-consistent mode-coupling theory for fore we argue that 0 r p-p and focus on the cal-
the viscosity of unentangled polyelectrolyte solutions. Start- culation of the polymer contribution. r is calculated by
ing with the polymer center-of-mass density as the slow vari- neglecting the effect of the interactions on the rotational de-
able, we develop an expression for the polymer contribution gree of freedom and is given by22
to the viscosity that is identical to that of Geszti.16 The in- 2 cp
termediate scattering function is obtained from the self- r r, 2.2
¨ 15 N
consistent mode coupling theory as in the approach of Gotze
and co-workers.17 The center-of-mass structure factor, re- where r is the rotational friction coefficient and is evaluated
quired in the theory, is obtained approximately from integral by a hydrodynamic calculation of an ellipsoid of the aspect
equation theory.18 –20 ratio N as23
The theory explains the behavior of the viscosity based 3
purely on the behavior of the static structure factor. For short 2 N4 1
r . 2.3
chains, the theory predicts that in salt-free solutions R de- 3 2N 2 1 2
creases monotonically with increasing concentration. With ln N N 1 N
N2 1
small amounts of added salt R displays a shallow peak as a
function of concentration at low concentrations. For longer For p-p , we employ a mode coupling treatment similar to
chains, the theory predicts a peak in R as a function of that of Geszti16 to derive a microscopic expression.
concentration for all salt concentrations, including salt-free The first step in the mode coupling approach is the
solutions. The peak occurs at concentrations slightly lower choice of the slow collective variables for the description of
than the overlap threshold concentration c * . With the addi- the dynamics of the required correlation functions. Natural
tion of salt, the intensity of the peak diminishes, but the choice is the hydrodynamic variables, i.e., the three momen-
position is unchanged. For a given salt concentration, the tum current densities of polyion, J (k), for the co-ordinates
height of the peak increases dramatically as the degree of x, y, and z, and the polyion number density P (k) de-
NP
polymerization is increased, but the position is unchanged. fined as P (k) i 1
e ik"ri where ri is the position of the
These predictions can be tested experimentally. center of mass of the ith polyion and N P is the number of
polyions. A calculable microscopic relation for the viscosity culations using stick boundary conditions. For the ellipsoid
is obtained by using the projection-operator formalism to with the aspect ratio of N, the total diffusion coefficient is
rewrite the well-known Green-Kubo time correlation func- given by23
tion expression in terms of P and Q operators.21 The stan-
dard approximation in the mode-mode coupling expansion is D 2D k BT 1
D0 ln N N2 1 . 2.8
to consider the subspace of various binary products of the 3 3 0 N2 1
basic slow variables. Among such binary products, the odd
Note that this theory considers only the translational mo-
ones with respect to time inversion do not contribute to the
tion of the rods, which is assumed to be isotropic. We argue
viscosity, and only the even combinations can be retained.
that we can neglect the anisotropy in translation and its cou-
The two obvious choices of the binary product are the
pling to rotation in the concentration regimes we consider.
density-density term and the current-current term. The cur-
We estimate the contribution from the anisotropy and cou-
rent terms are expected to decay much faster than the density
pling to diffusion as follows: In the dilute limit, if the inter-
term, due to the friction with the surrounding solvent mol-
action between polyions is neglected, it is possible to solve
ecules. Thus, we neglect this contribution. Finally, all four-
the rotation-coupled diffusion equation and evaluate F(k,t)
particle correlations are approximated as the product of two-
exactly.22,26 –29 If the ratio between the parallel (D ) and per-
particle correlations. With the above approximations and
pendicular (D ) diffusion coefficients is not very large, the
simplifications, the final expression for the zero frequency
change in the relaxation rate of F(k,t) at short times arising
viscosity is written in terms of the density correlation func-
from a coupling with rotational diffusion is also small. For
tion of the polyion as
example, for D /D 2, the relaxation rate of F(k,t) at
k BT S k 2
F k,t 2 short times is changed by less than 10% due to the coupling
0 r dkk 4 dt , with the rotational diffusion. At longer times,
60 2 0 0 S k S k
t R , where R is the rotational relaxation time, F(k,t) is
2.4
simply given by exp D0k2t , where D 0 is the average dif-
where S(k) is the static center-of-mass structure factor of fusion coefficient defined by D 0 (D 2D )/3. Thus, de-
the polyions, F(k,t) is the corresponding intermediate scat- coupling rotation from translation and assuming the transla-
tering function, S (k) is the derivative of S(k) with respect tion diffusion is isotropic are reasonable approximations in
to k. dilute solutions. These approximations become questionable
In order to evaluate the viscosity, we need the interme- when the concentration or rod length becomes large when,
diate scattering function F(k,t) for the polyions and it due to the entanglement effects, the rotational time increases
should also be evaluated using MCT. As discussed above, we steeply and anisotropy of the translational diffusion will be
again assume that dynamics of counterions and solvent mol- enhanced. This regime, however, is far beyond the scope of
ecules is decoupled from dynamics of polyion. Then, the the present paper.30
equation for F(k,t) is expressed in a closed form as17
III. STATIC PROPERTIES
D 0k 2 t
F k,t F k,t dt M k,t t To proceed further we require a model for the polyions
t S k 0
and a means of calculating the static structure of the polyion
centers of mass. In this work the molecules are modeled as a
F k,t , 2.5 collection of interaction sites arranged linearly in a rodlike
t
configuration. Each particle consists of N tangent charged
where D 0 is the bare collective diffusion coefficient. hard spheres or sites with hard sphere diameter , which is
M (k,t) is the memory kernel given by used as the unit of length in this paper. Each sphere carries a
negative fractional charge f e, where e is the charge on an
PD 0 dq 2 electron. The effect of solvent and small ions is included into
M k,t V q,k q F k q ,t F q,t ,
2 2 3 k the potential of interaction between sites on the polyelectro-
2.6 lyte molecules. The resulting effective potential u(r) is
given by
where P P (k 0) c p /N is the average number den-
sity of polyion and V(q,k q) is the vertex function given in for r
terms of the direct correlation function c(q) as u r 3.1
exp r /r for r ,
V q,k q k•qc q
ˆ k• k q c k q .
ˆ 2.7 where 1/k BT, f 2 l B /(1 ), l B e 2 / is the
Bjerrum length, is the dielectric constant of the solvent,
Equation 2.5 is a standard MCT equation familiar in the and is the inverse screening length,
supercooled liquids and colloids17,24,25 community. This 4 l B ( f 2 c p 2c s ) where c s is the number density of the
equation is a nonlinear integro-differential equation which monovalent salt, and c p is the number density of polymer
has to be solved self-consistently. The numerical procedure sites. In all the calculations presented in this work, l B
to solve Eq. 2.5 is elucidated in Ref. 24. 0.758 and f 1. If 4 Å, then an added salt concen-
The bare diffusion coefficient D 0 is obtained from the tration of 1 mM corresponds to a reduced salt concentration
value for a long ellipsoid calculated from hydrodynamic cal- of c s 3 4 10 5 .
FIG. 1. Comparison of PRISM predictions solid lines for the center-of-
mass structure factor to Monte Carlo simulation results for c p 3 10 3 ,
lB , and N 20. Dashed lines are simulation results for S ss (k)/ ˆ (k). FIG. 2. Static structure factor S(k) predicted using the PRISM theory solid
lines and the DH approximation dotted lines for N 150 and c s
1 mM, and for various polyion concentrations. From left to right, the
polyion concentrations are c p 3 10 6 , 2 10 5 , 10 4 , and 5 10 4 .
The center-of-mass static structure factor is calculated
using integral equations. The single chain structure factor
ˆ (k) is known exactly for this model. The site-site static
tion for various concentrations of added salt, and for degrees
structure factor S ss (k) is obtained from the polymer refer-
of polymerization of N 20 and 150, respectively. The
ence interaction site model PRISM theory,31 as described
added salt concentrations of 1, 2, and 5 mM correspond to
elsewhere.19 It has previously been established,32 by direct
reduced salt concentrations of c s 3 4 10 5 , 8 10 5 ,
comparison of theoretical predictions for S ss (k) to computer
and 2 10 4 , respectively. In the figures, the abscissa is the
simulations that PRISM is accurate for S ss (k). The center-
polymer concentration divided by the overlap threshold con-
of-mass structure factor is the approximated as S(k)
centration c * which, for this model, is given by c * 3
S ss (k)/ (k). To check the validity of this approximation,
1/N 2 . In all cases we find that the major contribution
we perform Monte Carlo simulations of rods interacting via
comes from the polymer-polymer interaction given by the
screened Coulomb interactions, and calculate S ss (k) and
mode-coupling expression in Eq. 2.4 . The contribution of
S(k). The simulation algorithm is identical to that described
the rotational Brownian motion of the individual rod r is
elsewhere32 except that we do not perform the Ewald sum.
independent of the polyion density and does not affect the
Figure 1 compares simulations results for S(k) filled circles
qualitative behavior except for the low concentration regime
and S ss (k)/ (k) dotted lines for l B , c p 3 10 3 , and
where the mode-coupling contribution becomes very small.
N 20, and shows that the approximation for S(k) is quite
For N 20 Fig. 3 a and c s 0 salt-free , R is a
accurate. Also shown in the figure is the PRISM prediction
monotonically decreasing function of c p /c * . As the salt con-
for S(k). The PRISM S(k) correctly reproduces the liquid-
centration is increased, the value of R decreases at all poly-
like structure manifested in the peak of S(k). In fact, the
mer concentrations. For N 20 this results in a shallow peak
theory is in quantitative agreement with the simulation re-
in R at low polymer concentrations and 1 mM salt. For high
sults.
values of added salt R is a monotonically increasing func-
Figure 2 depicts S(k) from PRISM for N 150 and for
tion of polymer concentration. These results are typical of
several polyion concentrations. We also show the results de-
cases when the static structure factor does not display a very
rived from a simpler Debye-Huckel DH approximation.
strong peak at low wave vectors. The influence of the long-
The DH result is derived by taking the →0 limit, and ap-
range liquidlike order on the dynamic properties is therefore
proximating the site-site direct correlation function, c ss (r),
very weak. The predictions for short chains are qualitatively
by c ss (r) u(r) for all r. The resulting intermolecular
similar to other theories6,9 that ignore the effect of static
structure factor, denoted S DH(k), is given by
structure on the dynamic properties.
1 As the degree of polymerization N is increased, the
S DH k . 3.2 theory predicts a prominent peak in R that occurs at a con-
4 lB
1 2 cp k centration just below the overlap threshold concentration.
k2
The amplitude of this peak increases with increasing degree
of polymerization and decreases with increasing salt concen-
IV. RESULTS
tration. This can be seen in Fig. 3 b which depicts R as a
For salt-free solutions, the theory predicts that the re- function of c p /c * for N 150. In salt-free solutions, the
duced viscosity R is a monotonically decreasing function of peak in R is very prominent. The addition of salt dramati-
polymer concentration, for short chains. As the chain length cally reduces the height of the peak, although a peak is
is increased, a peak in R is predicted, at concentrations clearly present even for high 5 mM salt concentrations.
slightly below the overlap threshold concentration. Figures Note that both axes are plotted on a logarithmic scale in
3 a and 3 b depict R as a function of polymer concentra- order to fit all the curves on the same figure.
 FIG. 4. Dependence of the reduced viscosity R on polymer concentration
FIG. 3. Dependence of the reduced viscosity R on polymer concentration for various degrees of polymerization, N 20 , 50 ( ), 100 * , and 150
for various salt concentrations, c s 0 mM , 1 mM ( ), 2 mM ( * ), and , and for a c s 1 mM and b c s 2 mM. Note that the abscissa is
5 mM , and for a N 20 and b N 150. Note that the abscissa is concentration divided by the overlap threshold concentration c * and both
concentration divided by the overlap threshold concentration c * , and both axes are logarithmic.
axes are logarithmic.
V. DISCUSSION
The reduced viscosity is a strong function of chain The main ingredients of the theory of this work are i
length, in a manner that depends on the salt concentration. the use of a fully self-consistent mode-coupling approxima-
Figures 4 a and 4 b depict R as a function of polymer tion SCMCT , and ii accurate estimates of the structure of
concentration for various values of N and c s 1 and 2 mM, the solution. It is of interest to determine how the actual
respectively. In both salt-free and added salt solutions, the predictions depend on these two components. We compare
peak in R grows with increasing degree of polymerization,
but the position of the peak is insensitive to the value of N.
Figure 5 depicts the value of R at the maximum as a func-
tion of degree of polymerization for various salt concentra-
tions. The peak values is fitted well by a power law except
for the salt-free case. The molecular weight dependence is
very strong, much stronger than what is obtained for en-
tangled neutral polymer melts. The exponent decreases dra-
matically as the salt concentration is increased.
The physical interpretation of these results is that the
peak in R arises from intermolecular correlations between
the rods. The main physical feature that is input into the
theory is an accurate estimate of the static structure factor of
the polyelectrolyte solutions. The viscosity is then calculated
using a fully self-consistent mode-coupling theory. Any scal-
ing analysis of the dependence of R on N and c p must take FIG. 5. N dependence of the peak value of R for various salt concentra-
into account the complex dependence of static correlations tions. Dotted lines are power law fits N , with 5.8, 3, 2.2, and 1.5, for
on the dynamics. c s 0, 1, 2, and 5 mM, respectively.
 B. Influence of static structure
In order to see how the accurate estimates of the struc-
ture affects the results, we compare MCT predictions for the
viscosity using PRISM results for the structure factor to
those using the DH approximation for S(k) given by Eq.
3.2 . A good test of the importance of liquid structure would
be to compare SCMCT with PRISM S(k) to SCMCT with
the DH S(k). We find, however, that SCMCT with the DH
S(k) predicts a so-called ergodic-nonergodic transition, i.e.,
F(k,t) fails to relax to zero, for low polyion concentrations.
Such a transition, which is also predicted by the MCT for
supercooled liquids, leads to a divergence in the viscosity at
FIG. 6. Comparison of predictions for the concentration dependence of R
finite concentrations! This prediction is clearly incorrect,
from a fully self-consistent MCT SCMCT — and lowest order MCT since no such divergence is seen in experiment, or expected
LMCT for N 150. PRISM results for S(k) are used in all on physical grounds. This emphasizes, however, that
calculations. In each case, the curves correspond to from top to bottom salt SCMCT is very sensitive to the structure factor used as in-
concentrations of c s 0 mM , 1 mM ( ), 2 mM * , and 5 mM .
put, as one would expect. We attribute the fictitious transition
to an overestimation of the memory kernel of Eq. 2.6 at
the predictions of this work with those of related, and sim- large wave vectors which allows anomalous positive feed-
pler, theoretical schemes, as described below. back into the relaxation of F(k,t). This over estimation of
the memory kernel arises from the broadened and featureless
A. Comparison of SCMCT with lowest order MCT S(k) in the DH approximation. In reality, the hard-core in-
„LMCT… teraction comes into play at large k, thus resulting in a flat-
tening of S(k). The hard-core interaction is, of course, ne-
The MCT requires an expression for the intermediate glected in the DH approximation.
scattering function F(k,t), which we obtain from the self- We therefore investigate the influence of structure using
consistent mode coupling equation, Eq. 2.5 . A simpler ap- the LMCT. A combination of the DH and LMCT approxima-
proximation amounts to neglecting the memory kernel tions allows us to derive simple scaling results for the vis-
M (k,t) in Eq. 2.5 . We refer to this approximation, where cosity. Since the form factor (k) depends on k only weakly,
F(k,t) is simply given by we can set (k) (0) N in Eq. 3.2 to get
D 0k 2t
F k,t S k exp 5.1 1
S k S k 3
4 l B Nc p
as the lowest order MCT LMCT . Such an approximation 1
k2 2
has been previously investigated by others6,3 but with the
¨
Debye-Huckel DH approximation for the static structure. 1
, 5.3
The theory of Nishida et al.9 is closely related to the LMCT 1 A x
with an approximate concentration independent structure
where x k/ , A 4 l B Nc p 3 / 2 Nc p 3 /(c p 2c s ), and
factor obtained from numerical calculations at zero density.
(x) 1/(1 x 2 ). With these simplifications, R from the
In the LMCT, Eq. 2.4 is readily integrated over time to give
LMCT is given by
2
k BT S k
R,LMCT 2 dkk 2 . 5.2 k B TA 2 x 2
120 0c pD 0 0 S k R dxx 2 2, 5.4
120 2 0 c p 3