Statics, Dynamics,and Rheological Propertiesof Micellar ...

pawunderarmΜηχανική

18 Ιουλ 2012 (πριν από 5 χρόνια και 1 μήνα)

661 εμφανίσεις

Statics, Dynamics, and Rheological Properties of
Micellar Solutions by Computer Simulation
Thèse en cotutelle avec l’
U
niversité
L
ibre de
B
ruxelles
THÈSE
présentée à
l’Université Paul Verlaine-Metz
par
Chien-Cheng HUANG
pour l’obtention du grade de Docteur
Spécialité : Physique
Soutenue le 25 septembre 2007

devant la commission d’examen
M. Wim

Briels
Professeur
,
Twente University, Pays-Bas
Rapporteur
M. Roland G. Winkler
Professeur, Forschungszentrum Jülich, Allemagne Rapporteur
M. Michel Mareschal
Professeur
,
Université Libre de Bruxelles
Examinateur
M. Joachim Wittmer
D
R CNRS
,
Institut Charles Sadron, Strasbourg
Examinateur

Mme

Hong Xu
Professeur
,
Université Paul Verlaine-Metz
Directeur de Thèse
M. Jean-Paul Ryckaert
Professeur
,
Université Libre de Bruxelles
Directeur de Thèse
Remerciements

Dans le cadre d’une thèse en co-tutelle entre le Laboratoire de Physique des Milieux
Denses de l’Université Paul Verlaine-Metz et le Laboratoire de Physique des
Polymères de l’Université de Libre de Bruxelles, j’ai été amené à travailler au sein de
deux laboratoires, ce qui fut une expérience très enrichissante tant d’un point de vue
scientifique qu’humain.
J’exprime, tout d’abord, ma reconnaissance envers le Prof. Hong Xu pour sa
pédagogie, sa grande disponibilité et son soutien ainsi que pour la formation qu’elle
m’a dispensé. Je tiens à remercier ensuite le Prof. Jean-Paul Ryckaert pour toutes les
discussions productives que nous avons eues ensemble et le partage de son savoir
dans le domaine de la physique. Je les remercie pour leur soutien jour après jour et
pour leur qualité humaine.

Merci à vous deux de m’avoir donné un si bel exemple de rigueur et d’honnêteté
scientifique.
Je remercie également le Dr. Joachim Wittme et le Prof. Marc Baus pour les
nombreuses idées sur les voies à explorer ainsi que pour les discussions très
enrichissantes.
Je remercie grandement Monsieur Georges Destrée pour les aides sur les techniques
de simulations.
Je souhaite aussi remercier Monsieur François Crevel pour les nombreuses
discussions.
Je voudrais également remercier l’ensemble des membres du jury de thèse, et en
particulier le Prof. Roland G. Winkler et le Prof. Wim Briels pour avoir accepté la
charge de rapporteurs de cette thèse. Merci au Prof. Michel Mareschal d’avoir bien
voulu présider le jury et Dr. Joachim Wittme pour l’intérêt qu’il a porté à ce travail.

Enfin, je tiens également à remercier ma famille et mes proches pour leur confiance et
leur soutien constant.
Abstract
Statics,dynamics,rheology and scission recombination kinetics of self as-
sembling linear micelles are investigated at equlibriumstate and under shear
flowby computer simulations using a newly proposed mesoscopic model.We
model the micelles as linear sequences of Brownian beads whose space-time
evolution is governed by Langevin dynamics.A Monte Carlo algorithm
controls the opening of a bond or the chain-end fusion.A kinetic param-
eter ω modelling the effect of a potential barrier along a kinetic path,is
introduced in our model.For equilibrium state we focus on the analysis
of short and long time behaviors of the scission and recombination mech-
anisms.Our results show that at time scales larger than the life time of
the average chain length,the kinetics is in agreement with the mean-field
kinetics model of Cates.By studying macroscopic relaxation phenomena
such as the average micelle length evolution after a T-jump,the monomer
diffusion,and the zero shear relaxation function,we confirm that the ef-
fective kinetic constants found are indeed the relevant parameters when
macroscopic relaxation is coupled to the kinetics of micelles.For the non-
equilibrium situation,we study the coupled effects of the shear flow and the
scission-recombination kinetics,on the structural and rheological properties
of this micellar system.Our study is performed in semi-dilute and dynam-
ically unentangled regime conditions.The explored parameter ω range is
chosen in order for the life time of the average size chain to remain shorter
than its intrinsic (Rouse) longest relaxation time.Central to our analysis
is the concept of dynamical unit of size Λ,the chain fragment for which the
life time τ
Λ
and the Rouse time are equal.Shear thinning,chain gyration
tensor anisotropy,chain orientation and bond stretching are found to de-
pend upon the reduced shear rate β
Λ
= ˙γτ
Λ
while the average micelle size
is found to decrease with increasing shear rate,independently of the height
of the barrier of the scission-recombination process.
Contents
Introduction 1
1 Theoretical Framework 7
1.1 Statistical mechanics derivation of the distribution of chain lengths...7
1.2 A kinetic model for scissions and recombinations.............13
1.2.1 Macroscopic thermodynamic relaxation..............16
1.3 Linear viscoelasticity of unentangled micelles...............17
1.3.1 The Rouse model...........................18
1.3.2 The Rouse model implications on the viscous response of a monodis-
perse polymer solution........................19
1.3.3 The theory of Faivre and Gardissat and the viscoelasticity of
micelles................................20
2 The mesoscopic model of worm-like micelles 23
2.1 The potential [19]...............................24
2.2 The Langevin Dynamics...........................25
2.2.1 Nonequilibrium LD technique....................28
2.3 Monte-Carlo procedure............................29
3 Equilibrium Properties 35
3.1 Static properties...............................36
3.1.1 List of simulation experiments and chain length distribution...36
3.1.2 Chain length conformational analysis................38
3.1.3 Pair correlation function of chain ends and the distribution of
bond length..............................40
3.2 Dynamic properties..............................42
i
3.2.1 Kinetics Analysis...........................43
3.2.2 Analysis of the Monomer Diffusion.................65
3.2.3 Zero-shear stress time autocorrelation functions..........71
3.2.4 Macroscopic relaxation behavior..................73
3.3 General comments..............................74
4 Non-equilibrium properties 77
4.1 Collective rheological behavior.......................78
4.1.1 Orientation of the chains......................78
4.1.2 Viscosity...............................82
4.2 Chain length distribution and chain size dependent properties......85
4.2.1 Effect of shear flow on distribution of chain lengths........85
4.2.2 Saturation effect on orientational properties............87
4.3 Scission-recombination kinetics under flow.................92
5 Conclusions 99
Bibliography 105
ii
Introduction
In the field of self-assembling structures,supramolecular polymers are attracting
nowadays much attention [1,56].The terminology of supramolecular polymer applies
to any polymer-like flexible and cylindrical superstructure obtained by the reversible
linear aggregation of one or more type of molecules in solution or in melt.These
supramolecular polymers are typical soft matter systems and the chain length dis-
tribution which determines their properties,is very sensitive to external conditions
(temperature,concentration,external fields,salt contents,etc...).Wormlike micelles
are one of the most common type of supramolecular polymer [2].Self-assembled stacks
of discotic molecules [20] and chains of bifunctional molecules [21] are other examples.
All these examples differ by the nature of the intermolecular forces involved in the self-
assembling of the basic units,but they lead to a similar physical situation bearing much
analogy with a traditional system of polydisperse flexible polymers when their length
becomes sufficiently large with respect to their persistence length.The specificity and
originality of these supramolecular polymers comes from the fact that these chains are
continuously subject to scissions at random places along their contour and subject to
end to end recombinations,leading to a dynamical equilibrium between different chain
lengths species.
It is well known and schematically illustrated in figure 1 that some surfactant
molecules in solution can self-assemble and form wormlike micelles [2].As their mass
distribution is in thermal equilibrium,they are therefore sometimes termed “equilib-
rium polymers” (EP) [1].Similar system of equilibrium polymer are formed by liquid
sulfur [33,34,35],selenium [36] and some protein filaments [37]
The micellar solutions exhibit fascinating rheological behavior which has been re-
cently reviewed and discussed theoretically [56].Quite commonly,under shear flow
in a Couette cell,an originally isotropic micellar solution undergoes a shear-banding
1
transition [3],producing a two phases system spatially organized in a concentric man-
ner with both phases lying either close to the inner or to the outer cylinder of the
rheometer.As the viscosity of both phases is different,a velocity profile with two
slopes is observed through the gap.This non-equilibriumphase separation is the object
of numerous studies [56].The shear thinning and orientational ordering of wormlike
micelles [55] resembles the usual phenomenon observed in polymer solutions but,in
many micellar solutions close to the overlap concentration[4],one can also observe a
shear-thickening behavior whose microscopic origin is still a matter of debate nowa-
days.Much better understood is the trend of entangled micellar solutions to display a
Maxwell fluid rheological behavior [2].The simple exponential relaxation behavior has
been theoretically explained by the reptation-reaction model[56],taking into account
the scission-recombination mechanism by which local entanglements can be released by
chain scission.The theoretical treatment of the scission-recombination process within
the rheological theoretical approaches,is usually based on a simple mean-field (MF)
approach in which correlations between successive kinetic events are fully neglected.
The way to take into consideration such correlated transitions has been detailed by
O’Shaughnessy and Yu [12] in order to explain some high frequency features of the rhe-
ological behavior of micellar systems.As the latter kinetics model interprets the stress
relaxation as the result of local effects between successive correlated transitions where
chain ends produced by a scission recombine with each other after a small diffusive ex-
cursion of the chain ends,such a rheological behavior was called “diffusion-controlled”
(DC) by opposition to the standard MF model.
Computer simulations at the mesoscopic scale offers an interesting route towards
the understanding of the generic properties of wormlike micelles solutions on the basis
of a well controlled microscopic (or mesoscopic) model.The first accent within the
simulation approach was put on testing the chain length distribution for various algo-
rithms producing,on a lattice or in continuum space,temporary linear self assembling
structures.Earlier Monte Carlo simulations using an asymmetric Potts model were
performed to study static properties of equilibrium polymers [38,39,40].Rouault and
Milchev proposed a dynamical Monte Carlo algorithm [41],based on the highly efficient
bond fluctuation model (BFM) [42,43].In great detail,Wittmer and co-workers [13,22]
investigate the static properties of EP in dilute and semi-dilute solutions using the
same BFM.They confirmed scaling predictions of the mean chain length in dilute and
2
semi-dilute limits.Rouault [15] verified the dependence of mean chain length with con-
centration using an off-lattice Brownian dynamics simulation.Dynamics at equilibrium
and direct simulation of systems in shear flow were also much investigated over the last
ten years.Kroger and co-workers study EP at equilibrium [58] and under shear flow
[17] by Molecular Dynamics using the particular pair potential FENE-C model.This
model is a variant of the traditional pair potential of the LJ+FENE type used to model
(permanent) polymers,a pair potential which diverges both as r goes to zero (repulsive
forces) and at a distance r = R
c
where the corresponding attractive force between
neighboring monomers also diverges.The FENE-C potential is a truncated form of the
usual form which is set to a constant beyond a certain distance r = r
m
< R
c
value,
creating a finite potential well for a bounded pair with respect to the unbounded pair
(r > r
m
).In this way,the covalent bond can break or recombine as the pair distance
crosses the r
m
value.Using the same model,Padding and Boek [16] investigated the
recombination kinetics and the stress relaxation of wormlike micellar systems.They
found that at high concentrations,the kinetics is close to a diffusion-controlled mecha-
nism.Milchev and co-works [26] study micelle conformations and their size distribution
by an off-lattice microscopic model,to study solutions of EP in a lamellar shear flow
while Padding and Boek investigate the influence of shear flow on the formation of
rings in a EP system using FENE-C model [44].All these studies predict a decrease of
the average micelle size as a result of the shear flow.Other simulation studies envisage
that ultimately,various simulations at different length and time scale will have to be
combined to fill the huge gap between the atomistic length and time scales where the
precise chemistry is relevant and the mesoscopic scale where rheological properties can
really be probed by simulation.Studies of this type mixing both the atomistic and
mesoscopic approaches to study wormlike micelles rheology are under progress[57].
The aim of our thesis on isotropic wormlike micelles solutions is to analyze,on the
basis of a dynamical simulation at the mesoscopic level,the influence on the macro-
scopic relaxation phenomena of the dynamical coupling between the usual “flexible
polymer” relaxation processes and the “scission-recombination” kinetics.We will re-
strict ourselves to unentangled solutions (working slightly below or above the isotropic
semi-dilute threshold with a very flexible mesoscopic equilibrium polymer model) to
avoid the prohibitive computer time needed to follow relaxation phenomena governed
by entangled dynamics.The relevance of the MF kinetics model and the microscopic
3
origin of the scission-recombination rate constants is one important target of our study.
Unentangled supramolecular polymers dynamics was indeed found to be relevant in
selenium rheology.At the occasion of this experimental rheological study,Faivre and
Gardissat [36] proposed an extension of the traditional Rouse theory to include the
scission kinetic effects on the Rouse modes relaxation.They predicted the way the
zero shear rate viscosity decreases as the scission and recombination kinetic constants
increase,introducing the concept of dynamical subunit whose size fixes the relaxation
time which governs the shear modulus relaxation.
The thesis is organized as follows:The theoretical aspects of wormlike micelles
statics (distribution of micelle sizes) and the dynamics (Cates MF kinetic model,Faivre
and Gardissat theory[36])are gathered in chapter I.
In chapter II,we detail the particular mesoscopic Langevin Dynamics model which
is adopted throughout our thesis.Given our aims,this model has the particularly
useful feature that the static properties (in particular the distribution of chain lengths)
and the dynamics of the scission-recombination can be tuned separately so that we
will often investigate how the relaxation at a unique state point (from a static point of
view) is modified by tuning the kinetic rate constants.
Chapter III reports the results of a series of simulations at equilibrium,performed
by Langevin Dynamics.Both a dilute and a semi-dilute state points are treated.We
check our results regarding the distribution of chain sizes with respect to theoretical
prediction and analyze in detail the microscopic origin of the relevant rate constants.
Chapter IV reports a systematic study of a single semi-dilute state point of micellar
solution under shear flow.Here two kinds of parameters,the shear rate and the scission-
recombination rate constants,are systematically varied.We discuss the evolution of
the average micelle size with shear rate and relate it to a different evolution of the two
types of rate constants (scission and recombination).We also discuss the nature of
the relevant reduced shear rate and extrapolate our viscosity data to zero shear rate,
allowing us to test the Faivre-Gardissat theory’s predictions [36].
Our general conclusions are gathered in the last chapter V.
4
Cylindrical micelles
Equilibrium polymer
COARSE−GRAINING
Surfactant
k
k
s
r
B
q
Chain end
ScissionRecombination
E
Figure 1:Some surfactant molecules in solution self-assemble and form long wormlike
micelles which continuously break and recombine.Their mass distribution is,hence,in
thermal equilibriumand they present an important example of the vast class of systems
termed “equilibrium polymers”[1].The free energy E of the (hemispherical) end cap of
these micelles has been estimated [2] to be of order of 10k
B
T.This energy penalty
(together with the monomer density) determines essentially the static properties and
fixes the ratio of the scission and recombination rates,k
s
and k
r
.Additionally,these
rates are influenced by the barrier height B which has been estimated to be similar to
the end cap energy.Both important energy scales have been sketched schematically as
a function of a generic reaction coordinate q (see chapter 8 of reference [23]).Following
closely the analytical description [2,22,13] these micellar systems are represented in
this study by coarse-grained effective potentials in terms of a standard bead-spring
model.The end cap free energy becomes now an energy penalty for scission events,
i.e.,the creation of two unsaturated chain ends.The dynamical barrier is taken into
account by means of an attempt frequency ω = exp(−B/k
B
T).If ω is large,suc-
cessive breakage and recombination events for a given chain can be assumed to be
uncorrelated and the recombination of a newly created chain ends will be of standard
mean-field type.On the other hand,the (return) probability that two newly created
chain ends recombine immediately must be particularly important at large ω.These
highly correlated “diffusion controlled” [12] recombination events do not contribute to
the effective macroscopic reaction rates which determine the dynamics of the system.
6
Chapter 1
Theoretical Framework
In this chapter we will briefly present the theoretical framework of cylindrical micelles
and some standard polymer theoretical aspects which are pertinent to the flexible
micelles system.In Section 1,we introduce the statistical mechanics derivation of the
distribution of chain lengths and of the corresponding average chain length.Then in
section 2,for the scission-recombination kinetics at equilibrium,we study the theoretical
formalismof equilibriumpolymers (EP) developed by Cates and co-workers [9,2],based
on a mean-field approximation.And the third section of this chapter will focus on the
linear viscoelasticity of EP.In this section,a framework for linear viscoelacity of dilute
polymer solutions,and intrinsic shear modulus by the Rouse model are presented.
Finally,we briefly review the theory of Faivre and Gardissat [36] where a modification
of the standard Rouse theory of linear viscosity of a polydisperse polymer system is
proposed.
1.1 Statistical mechanics derivation of the distribution of
chain lengths
To treat a system of wormlike micelles theoretically,it is convenient to work at the
mesoscopic scale using a model of linear flexible polymers made of L monomers of size
b linked together by a non permanent bonding scheme.Within the system,individual
chain lengths fluctuate by bond scission and by fusion of chain ends of two different
chains.Statistical mechanics can be employed to predict the equilibrium distribution
of chain lengths[2].In terms of the equilibriumchain number density c
0
(L),the average
chain length L
0
and the total monomer density φ are given by
7
L
0
=
P

L=0
Lc
0
(L)
P

L=0
c
0
(L)
(1.1)
φ =

X
L=0
Lc
0
(L) =
M
V
(1.2)
where M denotes the total number of monomers in the system and V is the volume.
Conceptually,we consider the Helmholtz free energy F(V,T;{N(L)},N
s
) of a mix-
ture of chain molecules of different length L in solvent where,in addition to the temper-
ature T,the volume V and the number of solvent molecules N
s
,the number of chains
of each specific length N(L) is fixed.Let F(V,T;M,N
s
) be the Helmholtz free energy
of a similar system where only the total number of solute monomers M is fixed.The
equilibrium chain length distribution c
0
(L) = N(L)/V will result from the set {N(L)}
which satisfies the condition
F(V,T;M,N
s
) = min
{N(L)}
"
F(V,T;{N(L)},N
s
) +

X
L=0
LN(L)
#
(1.3)
The parameter  is the Lagrange multiplier associated with the constraint that
the individual number of chains N(L) must keep fixed the total number of monomers
M =
P

0
LN(L).Minimization requires that the first derivative with respect to any
N(L

) variable (L’=1,2,...) is zero,giving
δF(V,T;{N(L)},N
s
)
δN(L

)
+L

= 0;L

= 1,2,...(1.4)
We expect the entropic part of the total free energy F(V,T;{N(L)},N
s
) to be the
sumof translational and chain internal configurational contributions which both depend
upon the way the M monomers are arranged into a particular chain size distribution.
For the translation part,the polydisperse system entropy is estimated as the ideal
mixture entropy S
id
S
id
(V,T;{N(L)},N
s
) = −k
B
X
L
N(L) ln(CN(L)) +S
solv
(1.5)
where C = b
3
/V is a dimensionless constant independent of L and where S
solv
is
the solvent contribution,independent of the N(L) distribution.The configurational
entropy of an individual chain with L monomers is written as S
1
(L) = k
B
lnΩ
L
,in
8
1.1 Statistical mechanics derivation of the distribution of chain lengths
terms of Ω
L
,the total number of configurations of the chain.Adding the configurational
contributions to S
id
as given by eq.(1.5),the total entropy becomes
S(V,T;{N(L)},N
s
) = −k
B
X
L
N(L) [ln(CN(L)) −lnΩ
L
] +S
solv
(1.6)
We now turn to the energy E(V,T;{N(L)},N
s
) of the same system.If E
1
(V,T;
L) represents the internal energy of a chain of L monomers and E
s
the energy of a
solvent molecule,the energy can be written as
E(V,T;{N(L)},N
s
) =
X
L
N(L)E
1
(V,T;L) +N
s
E
s
(V,T) (1.7)
The key contribution in E
1
is the chain end-cap energy E which corresponds to
the chain end energy penalty required to break a chain in two pieces.We will suppose
that E
1
(L) = E + L˜ǫ where ˜ǫ is an irrelevant energy per monomer as M˜ǫ,its total
contribution to the system energy,is independent of the chain length distribution.
The present approximation of the total free energy of the system is thus given by
incorporating in the general expression (1.3) the expressions (1.6) and (1.7),giving
βF(V,T;{N(L)},N
s
) =
X
L
N(L) [lnN(L) +lnC −lnΩ
L
+βE +β˜ǫL] (1.8)
where irrelevant constant solvent terms have been omitted as we only need the first
derivative of the free energy with respect to N(L),which now takes the form
δβF
δN(L)
= lnN(L) +lnC −lnΩ
L
+βE +β˜ǫL+1.(1.9)
With this expression,the minimization condition on N(L) becomes
lnN(L) +lnC −lnΩ
L
+βE +1 +

L = 0 (1.10)
where 

= (β +β˜ǫ).We note at this stage that the second derivative of βF(V,T;
{N(L)}) +

P
L
LN(L) with respect to N(L) and N(L

) variables gives the non neg-
ative result
δ
LL

N(L)
,indicating that the extremum is indeed a minimum.
Solving for N(L) in eq.(1.10),we get
N(L) = C

−1
exp−(

L+βE −lnΩ
L
) (1.11)
where C

= eC while,according to eq.(1.2)),

must be such that
X
L
Lexp−(

L+βE −lnΩ
L
) = MC

(1.12)
9
The equilibrium N(L) variables are also related to the equilibrium chain length
average L
0
(see eq.(1.1)),so that
X
L
exp−(

L+βE −lnΩ
L
) =
MC

L
0
(1.13)
To progress,we nowneed to specify the explicit L dependence of Ω
L
.The traditional
single chain theories of polymer physics provide universal expression of Ω
L
in terms of
the polymer size,the environment being simply taken into account through the solvent
quality and the swollen blob size in the semi-dilute (good solvent) case.
1.1.0.1 The case of mean field or ideal chains
The basic mean-field or ideal chain model for a L segments chain gives
Ω
id
L
=

C
1
z
L

(1.14)
where z is the single monomer partition function and C
1
a dimensionless constant.
Adapting eq.(1.11),one has
N(L) =
C
1
C

exp−(βE) exp(−”L) (1.15)
where ” = 

−lnz must,according to eq.(1.12),be such that
X
L
Lexp−(”L) =
1
”
2
=
MC

C
1
exp−(βE)
(1.16)
while eq.(1.13) takes the form
X
L
exp−(”L) =
1
”
=
MC

L
0
C
1
exp−(βE)
(1.17)
In eqs.(1.16) and (1.17),sums over L from 1 to ∞have been approximated by the
result of their continuous integral counterparts.
Combining eqs.(1.15),(1.16) and (1.17),one gets the final expression for the chain
number densities
c
0
(L) =
φ
L
2
0
exp(−
L
L
0
) (1.18)
with the average polymer length given by
L
0
= B
1/2
1
φ
1
2
exp

βE
2

.(1.19)
where B
1
= eb
3
/C
1
is a constant depending upon the monomer size b and the prefactor
in the number of ideal chain configurations in eq.(1.14).
10
1.1 Statistical mechanics derivation of the distribution of chain lengths
1.1.0.2 The case of dilute chains in good solvent
Polymer solutions are in a dilute regime when chains do not overlap and in semi-dilute
regime when chains do strongly overlap while the total monomer volume fraction is still
well below its melt value.In the semi-dilute regime in good solvent condition,chains
remain swollen locally over some correlation length,known as the swollen blob size χ,
but they are ideal over larger distances as a result of the screening of excluded volume
interactions between blobs.
Specifically,for a given monomer number density φ,the blob size is given by the
condition that the blob volume χ
3
= b
3
L
∗3ν
is equal to the total volume V divided by
the total number of blobs M/L

in the swollen blob.This gives estimates in terms of
the reduced number density φ

= b
3
φ,
L

= φ

(
1
1−3ν
)
(1.20)
χ = bφ

(
ν
1−3ν
)
(1.21)
where ν = 0.588 in present good solvent conditions[8].
In living polymers characterized by a monomer number density φ and some averaged
chain length L
0
,the semi-dilute conditions correspond to the case L
0
>> L

.We
discuss in this subsection the theory for the dilute case where L
0
<< L

.We will come
back to the semi-dilute case in the next subsection.
Self-avoiding walks statistics apply to dilute chains in good solvent,and we thus
adopt the number of configurations[6,7](See especially page 128 of the book of Grosberg
and Khokhlov [7])
Ω
EV
L
= C
1
L
(γ−1)
z
L
(1.22)
for a chain of size L,where γ is the (entropy related) universal exponent equal to 1.165,
z is the single monomer partition function and C
1
a dimensionless constant.
Incorporating expression (1.22) in eq.(1.11),one gets
N(L) =
V
B
1
exp−(βE)L
(γ−1)
exp−(”L) (1.23)
where B was introduced in eq.(1.19) and where ” = 

−lnz must be fixed by eq.(1.2)
X
L
L
γ
exp−(”L) = B
1
φexp(βE) (1.24)
11
while eq.(1.13) takes here the form
X
L
L
(γ−1)
exp−(”L) =
B
1
φ
L
0
exp(βE) (1.25)
If L is treated as a continuous variable,eqs (1.24) and (1.25) can be rewritten in
terms of the Euler Gamma function satisfying Γ(x) = xΓ(x −1) as
Z

0
L
γ
exp−(”L)dL =
Γ(γ +1)
”
(γ+1)
= B
1
φexp(βE) (1.26)
Z

0
L
(γ−1)
exp−(”L)dL =
Γ(γ)
”
γ
=
B
1
φ
L
0
exp(βE) (1.27)
From eqs (1.26) and (1.27),one gets
” =
γ
L
0
(1.28)
B
1
φexp(βE) =
Γ(γ +1)
γ
(γ+1)
L
(γ+1)
0
(1.29)
These results lead then finally to the Schulz-Zimmdistribution of chain lengths,namely
c
0
(L) =
exp(−βE)
B
1
L
(γ−1)
exp(−γ
L
L
0
) (1.30)
and an average polymer length given by
L
0
=

γ
γ
Γ(γ)

1
1+γ
B
1
1+γ
1
φ
1
1+γ
exp

βE
(1 +γ)

.(1.31)
1.1.0.3 The semi-dilute case
We consider here the semi-dilute case in good solvent where the average length of living
polymers L
0
is much larger than the blob length L

.The usual picture of a semi-dilute
polymer solution is an assembly of ideal chains made of blobs of size χ.Using this
approach,Cates and Candau [2] and later J.P.Wittmer et al [13] derived the relevant
equilibrium polymer size distribution.In this subsection,we adapt their derivation to
the theoretical framework presented above.
Let Ω
b
be the number of internal configurations per blob and z

some coordination
number for successive blobs.As there are n
b
= L/L

blobs for a chain of L monomers,
we write the total number of internal configurations of a chain of size L as
Ω
SD
L
= C
1
L
∗(γ−1)
Ω
L/L

b
z

L/L

(1.32)
12
1.2 A kinetic model for scissions and recombinations
where γ is the universal exponent in the excluded volume chain statistics met earlier
for chains in dilute solutions.The important factor L
∗(γ−1)
can be seen as an entropy
correction for chain ends just like E was an energy correction to L˜ǫ.This entropic
term which involves the number of monomers per blob,is needed to take into account
that when a chain breaks,its two ends are subject to a reduced excluded volume
repulsion.The other factors in eq.(1.32) will lead to terms linear in L after taking
the logarithm and thus will be absorbed in the Lagrange multiplier definition,as seen
earlier in similar cases for ideal and dilute chains.The resulting expression of N(L) in
terms of the Lagrange multiplier (cf.eq.(1.15)) can then be written by analogy as
N(L) =
C
1
C

exp−[βE −(γ −1) lnL

] exp(−L) (1.33)
Proceeding as in the ideal case (simply replacing at every step the constant βE
by βE − (γ − 1) lnL

,one recovers in the semi-dilute case the simple exponential
distribution
c
0
(L) =
φ
L
2
0
exp(−
L
L
0
) (1.34)
with a slightly different formula for the average polymer length
L
0
∝ φ
α
exp

βE
2

(1.35)
where α =
1
2
(1 +
γ−1
3ν−1
) is about 0.6.
1.2 A kinetic model for scissions and recombinations
The interest for wormlike micelles dynamics came from the experimental observation
that entangled flexible micelles often display,after an initial strain,a simple exponential
stress relaxation.The mechanism of such relaxation is different from that of usual
dead polymer entangled melts.In the latter system,stress relaxation requires that
individual chains leave the strained topological tube created by the entangled temporary
network by a reptation mechanism.A theoretical model,taking into account the extra
relaxation mechanismcaused by scissions and recombinations of micelles,leads [9] to an
exponential decay of the shearing forces with a decay Maxwell time equal to τ ≈

τ
b
τ
rep
where τ
rep
is the chain reptation time and τ
b
is the mean life time of a chain of average
size in the system.
13
We will assume in the following that the Cates’s scission-recombination model gov-
erning the population dynamics,originally devised to explain entangled equilibrium
polymer melt rheology,should also apply to the kinetically unentangled regime which
is explored in the present work.
This Cates kinetic model [9] assumes that
• the scission of a chain is a unimolecular process,which occurs with equal proba-
bility per unit time and per unit length on all chains.The rate of this reaction is
a constant k
s
for each chemical bond,giving
τ
b
=
1
k
s
L
0
(1.36)
for the lifetime of a chain of mean length L
0
before it breaks into two pieces.
• recombination is a bimolecular process,with a rate k
r
which is identical for all
chain ends,independently of the molecular weight of the two reacting species they
belong to.It is assumed that recombination takes place with a new partner with
respect to its previous dissociation as chain end spatial correlations are neglected
within the present mean field theory approach.It results from detailed balance
that the mean life time of a chain end is also equal to τ
b
.
Let c(t,L) be the number of chains per unit volume having a size L at time t.On
the basis of the model,the following kinetic equations can be written [9]
dc(t,L)
dt
= −k
s
Lc(t,L) +2k
s
Z

L
c(t,L

)dL

+
k
r
2
Z
L
0
c(t,L

)c(t,L −L

)dL

−k
r
c(t,L)
Z

0
c(t,L

)dL

(1.37)
where the two first terms deal with chain scission (respectively disappearance or appear-
ance of chains with length L) while the two latter terms deal with chain recombination
(respectively provoking the appearance or disappearance of chains of length L).
It is remarkable that the static solution of this empirical kinetic model leads to an
exponential distribution of chain lengths.Indeed,direct substitution of solution c
0
(L)
in the above equation leads to the detailed balance condition:
φ
k
r
2k
s
= L
2
0
(1.38)
14
1.2 A kinetic model for scissions and recombinations
the ratio of the two kinetic constant being thus restricted by the thermodynamic state.
Detailed balance means that for the equilibrium distribution c
0
(L),the number of
scissions is equal to the number of recombinations per unit time and volume.The total
number of scissions and recombinations per unit volume and per unit time,denoted
respectively as n
s
and n
r
,can be expressed as
n
s
= k
s
φ
L
2
0
Z

0
Lexp(−L/L
0
)dL = k
s
φ (1.39)
n
r
=
k
r
2
φ
2
L
4
0
Z

0
dL

Z

0
dL”exp(−
L

L
0
) exp(−
L”
L
0
) =
φ
2
k
r
2L
2
0
(1.40)
and it can be easily verified that detailed balance condition equation 1.38 implies n
s
=
n
r
.
Mean field theory assumes that a polymer of length L will break on average after
a time equal to τ
b
= (k
s
L)
−1
according to a Poisson process.This implies that the
distribution of first breaking times (equal to the survival times distribution) must be
of the form
Ψ(t) ∝ exp(−
t
τ
b
) (1.41)
for a chain of average size.Detailed balance then requires that the same distribution
represents the distribution of first recombination times for a chain end[2].Accordingly,
throughout the rest of this chapter,the symbol τ
b
will represent as well the average
time to break a polymer of average size or the average time between end chain recombi-
nations.In the same spirit,we stress that among the different estimates of τ
b
proposed
in this work,some are based on analyzing the scission statistics while others are based
on the recombination statistics.
Two additional points may be stressed at this stage:
• The mean field model in the present context has been questioned [12] because
in many applications,there are indications that a newly created chain end often
recombines after a short diffusive walk with its original partner.In that case,a
possibly large number of breaking events are just not effective and the kinetics
proceeds thus differently.
• Given the statistical mechanics analysis in the previous subsection,we see that
the equilibrium distribution of chain lengths resulting from the simple empirical
kinetic model is perfectly compatible with the equilibriumdistribution in polymer
15
solutions at the θ point (ideal chains) or for semi-solutions (ideal chains of swollen
blobs).The kinetic model is still pertinent,given its simplicity,in the case of
dilute solutions in good solvent as the chain length distribution although non
exponential,is not very far from it.
1.2.1 Macroscopic thermodynamic relaxation
In equation (1.35) the average length L
0
depends on the thermodynamic variables of
the systemsuch as the temperature,the pressure,and the volume fraction of monomers
φ.Under a sudden change in one of the thermodynamic variables (for example,the
temperature),the distribution of the length of the chains in equation (1.37) relaxes
to a new equilibrium exponential distribution characterized by a new average length
L

.The corresponding theory of macroscopic thermodynamic relaxation has been
given by Marques and Cates [27] and tested using a Monte Carlo study by Michev
and Rouault [46].This is in general a complicated,nonlinear decay which can be
monitored experimentally by light scattering which probes the evolution of the average
length.The theory is based on the equations of distribution of length of chain (1.34)
and the kinetic model of Cates [9] i.e.equation (1.37).The characteristic time provides
information about the kinetics of the system.
Marques and Cates [27] show that for any amplitude of the perturbation which
conserves the total volume fraction of monomers (for instance,a temperature modifica-
tion having an arbitrary time-dependence),the distribution function c(t,L) of equation
(1.37) remains exponential versus L.Indeed,they show that
c(t,L) = φf
2
(t) exp{−Lf(t)} (1.42)
is a non-linear eigenfunction of equation (1.37) with a eigenvalue f(t) that obeys
df
dt
= k
s

1 −
φk
r
2k
s
f
2
(t)

,(1.43)
where k
r
and k
s
the kinetic constants of the new equilibriumstate,after,e.g.a T-jump.
The time evolution of f(t) for a thermodynamic jump is then the solution of (1.43) with
the initial condition f(t = 0) = 1/L
0
.One obtains [27]:
f(t) = L
−1

tanh

t +t
0


if L
0
> L

(1.44)
f(t) = L
−1

coth

t +t
0


if L
0
< L

(1.45)
16
1.3 Linear viscoelasticity of unentangled micelles
where the constant t
0
is given by
t
0
= 2τ tanh

L

L
0

if L
0
> L

(1.46)
t
0
= 2τ coth

L

L
0

if L
0
< L

,(1.47)
and the time
τ =
1
2k
s
L

(1.48)
These equations describe the recovery of c(t,L) and hL(t)i the average chain length,
following a sudden change in temperature which is given by
hL(t)i =
1
f(t)
(1.49)
We see that by measuring the non-equilibrium τ the relaxation time of the mean chain
length,one can get an estimate of the microscopic scission rate constant k
s
.
1.3 Linear viscoelasticity of unentangled micelles
To describe the dynamics of cylindrical micelles in solution,each micelle can be mod-
elled as a linear collection of L identical fragments (L being proportional to the length
of the micelle),subject to Brownian motion to represent the dynamical effect of the
bath (the coupling to the rest of the degrees of freedom),subject to connecting forces
between adjacent fragments to enforce the linear structure of the micelle and possibly
subject to pair interactions between non connected fragment pairs.If,as the rela-
tive distance between adjacent fragments increases,the connecting tension force is not
bounded to a finite value to allow an intrinsic possibility of scission (e.g.the FENE-C
model of Kroger[17]),the scission and recombination processes must be incorporated
independently in the model.The latter must specify the way through which individ-
ual bonds can break (chain scission) or reform (chain end recombinations) by a pair
potential swap for the involved pair of fragments.The micelle dynamics thus needs
to combine all these various ingredients,implying that individual chain entities will
survive for a finite life time during which ordinary polymer relaxation takes place.
Within a theoretical approach of the dynamical properties of micellar solutions
where the emphasis lies in the search for analytical predictions,it is useful to explore
the validity of the traditional Rouse model to represent the “dead polymer” dynamics
17
part (valid as long as an entity survives),when it can be assumed that excluded volume
and entanglement effects are absent or weak.In the present purely theoretical discus-
sion,we thus restrict ourselves to the unentangled dynamics of θ solvent chains.We
briefly review the basic ingredients and essential results of the Rouse model for “dead
chains”,including linear viscoelasticity and we report the extension of this theory to
polymer chains when they are the object of scissions and recombinations,along the
lines proposed by Faivre and Gardissat [36].
1.3.1 The Rouse model
The Rouse model treats all chains as being independent.Consider a “dead chain” of
N elements at temperature T,where each fragment is subject to a friction force pro-
portional to its instantaneous velocity with friction coefficient ξ and to a random force
modelled as a 3D white noise stochastic process whose statistical properties are given
by the fluctuation-dissipation theorem in terms of T and ξ.The only interactions con-
sidered in the Rouse model are harmonic spring forces with spring constant k between
neighboring beads to ensure the linear connectivity.Under this dynamics,the average
squared distance between adjacent monomers defines a typical mean distance b given
by b
2
= 3k
B
T/k,where k
B
is the Boltzmann constant.The mean squared end-to-end
distance is simply (N −1)b
2
.
Chain dynamics considered via the end-to-end vector time relaxation function leads
at long times to an exponential decay with the main (Rouse) relaxation time
τ
R
= τ
0
N
2
(1.50)
where τ
0
=
ξb
2
3k
B

2
,indicating that the global relaxation time of a chain of size N is
proportional to N
2
.
The way to derive this important result requires the introduction of modes (Rouse
modes),defined as
¯
X
q
(t) =
1
N
N
X
j=1
A
qj
¯
R
j
(t);q = 0,1,2,   (N −1) (1.51)
where
¯
R
j
is the coordinate of the j
th
monomer in a chain and
¯
¯
A a square matrix with
elements
A
qj
= cos


N
(j −
1
2
)

;q = 0,1,2,   (N −1);j = 1,2,   N (1.52)
18
1.3 Linear viscoelasticity of unentangled micelles
According to this definition,the mode q = 0 represents the center of mass diffusive
motion,while the other (N-1) modes q > 0 are associated to the collective motion of
sub-chains of size N/q.The time autocorrelation of these internal chain modes are
<
¯
X
q
(t)
¯
X
q
(0) >=< X
2
q
> exp


t
τ
q

(1.53)
where
τ
q
=
ξb
2
3k
B

2

N
q

2
(1.54)
which shows that the Rouse time is the relaxation time associated to the mode q = 1,
while the relaxation time of higher index modes decreases as 1/q
2
.For completeness,
the mode amplitudes are given by
< X
2
q
>=
Nb
2

2
1
q
2
(1.55)
Rouse theory leads to an expression of the time autocorrelation of the chain end-
to-end vector
¯
S =
¯
R
N

¯
R
1
<
¯
S(t).
¯
S(0) >= (N −1)b
2
P′
q
τ
q
exp(−t/τ
q
)
P′
q
τ
q
(1.56)
where the prime in the sums means that the latter only involves the odd q mode
indices.This result shows a typical Rouse theory relaxation function:it is a sum
of terms implying the different mode relaxation,which explains that at long times,
only the slowest mode survives and the behavior of the time autocorrelation becomes
exponential with characteristic time τ
R
.
1.3.2 The Rouse model implications on the viscous response of a
monodisperse polymer solution
We now discuss the intrinsic shear modulus G(t) of a monodisperse polymer solution
which reflects the viscoelastic response of the system in the linear regime (limit of very
small shear rates).The Newtonian shear viscosity is related to the shear modulus
through
η
0
=
Z

0
G(t)dt (1.57)
Rouse theory for a monodisperse solution of “dead” polymers of size N at monomer
number density φ gives [48]
G(t) =
φ
N
k
B
T
N−1
X
q=1
exp(−2t/τ
q
) =
φ
N
k
B
T
N−1
X
q=1
exp(−2tq
2

R
) (1.58)
19
This important result indicates that,because of the independence of the dynamics
of the various chains in the system,Rouse theory gives for a collective and intensive
quantity like G(t) an expression which is proportional to the solute concentration while
its viscoelastic behavior is strictly identical to the single chain viscoelastic relaxation
process.Therefore,the shear modulus again appears as a sum of terms expressing the
specific contributions of the various single chain modes in the relaxation process,with
at long times,a simple exponential time decay with characteristic time τ
R
.By time
integration,according to eq.(1.57),the shear viscosity is then given (for N >> 1) by
η
0
=
π
2
12
φk
B
T
N
τ
R
(1.59)
where τ
R
is given in eq.(1.50).
1.3.3 The theory of Faivre and Gardissat and the viscoelasticity of
micelles
The above subsection is dealing with a monodisperse dilute solution of “dead” chains.
The stress relaxation function of a micellar system is in reality the object of a coupling
between the stress relaxation and the scission-recombination process if the time scale of
the second process is shorter than the viscoelastic times scales.In our work,we adapt
the theory proposed by Faivre and Gardissat [36],originally developed to interpret
rheological data of liquid selenium.Faivre and Gardissat [36] proposed a modification
of the standard Rouse theory of linear viscosity of a polydisperse polymer system [8]
to incorporate the influence of the scission events.
If we have a polydisperse system of polymers with normalized weight function W
p
,
the relaxation modulus given by the Rouse model should be given by
G(t) =

X
p=1
W
p
G
p
(t) (1.60)
where G
p
(t) concerns the “polymer” part of chain of size p which,within the Rouse
model,reads (see eq.(1.58))
G
p
(t) = G
0
1
p
p
X
q=1
exp
"
−2
t
τ
0

q
p

2
#
,(1.61)
where p is the polymerization degree,G
0
= φk
B
T a material constant,and τ
0
the
local dynamic relaxation time which depends on the solvent viscosity.Notice that for
20
1.3 Linear viscoelasticity of unentangled micelles
physical reasons we have q ≤ p in the sum of Rouse modes for any upper limit p.As
mentioned earlier,the stress relaxation is a superposition of exponentials,each term
corresponds to the contribution to the relaxation of a particular Rouse mode.All terms
have the same amplitude but decay with a characteristic time
τ
p
q
= τ
0
p
2
q
2
(1.62)
which depends on the p/q ratio which is the number of monomers per wave length for
that particular mode.If bond scission occurs independently and uniformly along the
chain with rate k
s
(per unit time and per bond),the lifetime of a chain of p/q monomers
should be
τ
b
=
q
k
s
p
.(1.63)
The Rouse mode q in the original chain of length p should therefore have a survival
probability which decays in time as exp(−k
s
(p/q)t),hence the idea of multiplying the
contribution of each mode to the relaxation of G(t) by this survival probability to take
into account the effect of bond scissions.This leads to the expression
G

p
(t) = G
0
1
p
p
X
q=1
exp
"
−k
s
p
q
t −2
t
τ
0

q
p

2
#
.(1.64)
When the life time τ
b
of the chain of average size is shorter than the internal Rouse
relaxation time of the same chain τ
0
L
2
0
,the viscoelastic response is independent of the
mean chain length L
0
but depends upon a new intrinsic time τ
Λ
which corresponds
to the relaxation of a dynamical unit of size Λ defined by the equality of the Rouse
relaxation time τ
0
Λ
2
and the life time,(k
s
Λ)
−1
:
τ
0
Λ
2
= (k
s
Λ)
−1
.(1.65)
It leads to
Λ = (τ
0
k
s
)
−1/3
(1.66)
and
τ
Λ
= τ
1/3
0
k
−2/3
s
(1.67)
so that equation (1.64) can be rewritten as
G

p
(t) = G
0
1
p
p
X
q=1
exp


t
τ
Λ
"
p

+2


p

2
#!
.(1.68)
21
Identifying the number of Rouse beads with the number of beads in present work,using
a known distribution of lengths of the chains of our system (1.18),the weight function
for our equilibrium polymers is
W
p
=
1
L
0
exp


p
L
0

.(1.69)
So that G(t) for the equilibrium polymers is now given by
G(t) =

X
p=1
W
p
G

p
(t).(1.70)
Note that this expression differs slightly from the Faivre and Gardissat final expres-
sion as apparently,in equation (18) of reference [36],they assume a relaxation time for
a segment of p monomers to be half of the usual Rouse time.(Note also that in their
paper,the symbol τ
0
corresponds to half the time τ
0
used in our work).
22
Chapter 2
The mesoscopic model of
worm-like micelles
The main goal our thesis is to exploit simulation techniques to improve our under-
standing of the link between the macroscopic behavior and microscopic aspects of the
structure and the dynamics of self-assembling micellar systems at equilibrium and un-
der shear flow.As it is pragmatically impossible to reach large enough scales of length
and time to determine those macroscopic properties using an atomistic level molecular
dynamics simulation approach,we adopt a mesoscopic model and a Langevin Dynamics
approach as done regularly in the studies on micellar systems.As sketched in figure 1,
micelles can be represented as linear sequences of Brownian beads which,in addition to
their usual Langevin Dynamics space-time evolution,can either fuse together to form
longer structures or break down into two pieces.The kinetics can be modelled by a
microscopic kinetic Monte-Carlo algorithm which generates new bonds between chain
ends adjacent in space or which breaks existing bonds between adjacent monomers.The
free energy E to creat this new end caps for these micelles becomes an energy penalty
for scission events,i.e.,the creation of two unsaturated chain ends.This scission energy
determines the static propeties and influences the scission and recombination rates,k
s
and k
r
.In this model,the scission energy E of micelles is controlled by a scission en-
ergy parameter W which is defined as an additive potential term which influences the
switching probability between bonded and unbounded potentials.The barrier energy
of recombination B also influences the kinetic rate constants.The advantage of our
model is that the dynamical barrier height B can be taken into account by means of
the attempt frequency ω associated to the Monte Carlo potential swap.If ω is small,
23
successive breakage and recombination events for a given chain can be assumed to be
uncorrelated and the recombination of newly created chain ends will be of standard
mean-field type.On the other hand,the (return) probability that two newly created
chain ends recombine immediately must be particulary important at large ω.These
highly correlated “diffusion controlled” [12] recombination events do not contribute do
the effective macroscopic reaction rates which determine the dynamics of the system.
The mesoscopic model is based on a standard polymer model with repulsive Lennard-
Jones interactions between all monomer pairs and an attractive FENE (finitely extensi-
ble nolinear elastic) potential [24] to enforce (linear) connectivity.For the scission and
recombination events,a Monte Carlo procedure is set up to switch back and forth be-
tween the bounded potential and the unbounded potential.This new model is thus at
the same time justified by its link with a standard polymer model and by the advantage
that the scission/recombination attempt frequency is a control parameter governing dy-
namics,without affecting the thermodynamic and structural properties of the system.
2.1 The potential [19]
We consider a set of micelles consisting of (non-cyclic) linear assemblies of Brownian
particles.Within such a linear assembly,the bond potential U
1
(r) acting between ad-
jacent particles is expressed as the sum of a repulsive Lennard-Jones (shifted and trun-
cated at its minimum) and an attractive part of the FENE type [24].The pair potential
U
2
(r) governing the interactions between any unbounded pair (both intramicellar and
intermicellar) is a pure repulsive potential corresponding to a simple Lennard-Jones
potential shifted and truncated at its minimum.This choice of effective interactions
between monomers implies good solvent conditions.
Using the Heaviside function Θ(x) = 0 or 1 for x < 0 or x ≥ 0 respectively,explicit
expressions (see figure 2.1) are
U
2
(r) = 4ǫ

(
σ
r
)
12
−(
σ
r
)
6
+
1
4

Θ(2
1/6
σ −r)) (2.1)
U
1
(r) = U
2
(r) −0.5kR
2
ln

1 −[
r
R
]
2

−U
min
−W (2.2)
In the second expression valid for r < R,k = 30ǫ/σ
2
is the spring constant and
R = 1.5σ is the value at which the FENE potential diverges.U
min
is the minimum
value of the sum of the two first terms of the second expression (occurring at r
min
=
24
2.2 The Langevin Dynamics
0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5
r
-10
0
10
20
30
40
50
60
70
80
U
Γ
W
0.96 < r < 1.2
Figure 2.1:Bounded potential U
1
(r) (continuous curve) and unbounded potential U
2
(r)
(dashed line) between a pair of monomers.W is a parameter tuning the energy required
to open the bond.The figure also shows the Γ region where potential swaps (corre-
sponding to bond scissions or bond recombinations) are allowed during the Brownian
dynamics simulation.The unit of length is σ,and for the energy is ǫ.
0.96094σ for the adopted parameters) while W is a key parameter which corresponds
to the typical energy gain (loss) when an unbounded (bounded) pair is the object of a
recombination (scission).
2.2 The Langevin Dynamics
The Langevin equation describes the Brownian motion of a set of interacting particles.
The action of the fluid is split into a slowly evolving viscous force and a rapidly fluc-
tuating random force.For a free particle in one dimension the equation is expressed
as
m
i
˙
¯v
i
(t) = −ξ¯v(t)
i
+
¯
F
i
(t) +
¯
R
i
(t) (2.3)
where barv
i
is the velocity of particle i with mass m
i
at time t,
¯
F
i
the systematic force,
ξ the friction coefficient,and
¯
R
i
represents the sum of the forces due to the incessant
collision of the fluid molecules.It is regarded as a stochastic force and is satisfying the
25
condition (fluctuation-dissipation theorem) [47]:
< R

(t)R

(t

) >= 2ξk
B

ij
δ
αβ
δ(t −t

) (2.4)
where the greek-letter subscripts αβ refer to the x,y or z components.
2.2.0.1 Brownian Dynamics (BD)
If we assume that the time scales for momentum relaxation and position relaxation
are well separated,then it is possible to consider only time intervals longer than the
momentumrelaxation times.In the diffusive regime the momentumof the particles has
relaxed and equation (2.3) can be modified by setting to zero the momentum variation
(left hand side of equation (2.4)).Then the “position Langevin equation” is given as
d¯r
i
dt
=
¯
F
i
(t)
ξ
+
¯
R
i
(t)
ξ
(2.5)
Earlier simulations within the present work were executed using an algorithm based on
(2.5) [25].The single BD step of particle i subject to a total force
F
i
is simply
¯r
i
(t +Δt) = ¯r
i
(t) +
¯
F
i
ξ
Δt +
¯
Δ
i
(Δt) (2.6)
where the last term given by
¯
Δ
i
(Δt) =
Z
t+Δt
t
¯
R
i
(t

)dt

/ξ (2.7)
corresponds to a vectorial random Gaussian quantity with first and second moments
given by
< Δ

(Δt) > = 0 (2.8)
< Δ
iαβ
(Δt)Δ

(Δt) > = 2(
k
B
T
ξ
)Δtδ
αβ
δ
ij
(2.9)
for arbitrary particles i and j and where αβ stand for the x,y,z Cartesian components.
At this stage,it is useful to fix units.In the following,we will adopt the size of
the monomer σ as unit of length,the ǫ parameter as energy unit and we will adopt
ξσ
2
/(3πǫ) as time unit.We also introduce the reduced temperature k
B
T/ǫ = T

,so
that in reduced coordinates (written here with symbol*) the random displacement Δ


is
< Δ


> = 0 (2.10)
< Δ


Δ


> =
2

Δt

δ
αβ
δ
ij
T

(2.11)
26
2.2 The Langevin Dynamics
which means that each particle,if isolated in the solvent,would diffuse with a RMSD of
p
2T

/π per unit of time.In the following,all quantities will be expressed in reduced
coordinates without the ∗ symbol.
2.2.0.2 Langevin Dynamics velocity Verlet scheme
Most of the simulation of the present work are performed using Langevin Dynamics
(LD) velocity Verlet algorithmes based on the work of [50].The explict expressions are
¯r
i
(t +Δt) = ¯r
i
(t) +c
1
¯v
i
(t)Δt +c
2
¯
F
i
(t)
m
Δt
2
+δ¯r
G
i
(2.12)
¯v
i(1/2)
= c
0
¯v
i
(t) +

c
0
c
2
c
1

¯
F
i
(t)
m
Δt +δ¯v
G
i
(2.13)
¯v
i
(t +Δt) = ¯v
i(1/2)
+
1
(γΔt)

1 −
c
0
c
1

¯
F
i
(t +Δt)
m
Δt (2.14)
where the first equation updates the position,the second equation updates the “half”
step velocities and the third equation completes the velocity move.γ = ξ/m,the
coefficients for the above equations are:
c
0
= e
−γΔt
(2.15)
c
1
= (γΔt)
−1
(1 −c
0
) (2.16)
c
2
= (γΔt)
−1
(1 −c
1
) (2.17)
c
3
= (γΔt)
−1

1
2
−c
1

(2.18)
where e
γΔt
is approximated by its power expansion for γΔt < 1.Each pair of vec-
torial component of δr
G
,δv
G
,i.e.δr
G

,δv
G

is sampled from a bivariate Gaussian
distribution[49] defined as
ρ

δr
G

,δv
G


=
1
2πσ
r
σ
v
(1 −c
2
rv
)
1/2
×exp
(

1
2(1 −c
2
rv
)


δr
G

σ
r

2
+

δv
G

σ
v

2
−2c
rv

δr
G

σ
r

δv
G

σ
v

!)
(2.19)
with zero mean values,variances given by
σ
2
r
=
D

δr
G


2
E
= Δt
2
k
B
T
m
(γΔt)
−1

2 −(γΔt)
−1

3 −4e
−γΔt
+e
−2γΔt


(2.20)
σ
2
v
=
D

δv
G


2
E
=
k
B
T
m

1 −e
−2γΔt

(2.21)
27
and the correlation coefficient c
rv
determined by
c
rv
=


δr
G

δv
G


= Δt
k
B
T
m
(γΔt)
−1

1 −e
−ξΔt

2
1
σ
r
σ
v
(2.22)
Each pair of cartesian components δr
G

and δv
G

are obtained by the appropriate Gaus-
sian distribution according to:
δr

= σ
r
η
1
(2.23)
δv

= σ
v

c
rv
η
1

2

p
1 −c
2
rv

(2.24)
where η
1
and η
2
are two independent random numbers with Gaussian distribution of
zero average and unit variance.
2.2.1 Nonequilibrium LD technique
In this sub-section we briefly present the technique adopted for imposing the shear flow
onto the system.We limit the study to a stationary,isovolumic and homogeneous planar
Couette flow characterized by a velocity field ¯u(r) =
¯
¯
k
T
r as illustrated in figure 2.2,
where the velocity gradient
¯
¯
k is constant in space and time.
Figure 2.2:Sketch of planar Couette flow.
For a planar Couette flow in the x direction with a shear along the y axis,the
velocity gradient k is defined as
¯
¯
k =






0 0 0
˙γ 0 0
0 0 0






where ˙γ is shear rate.With the imposition of the solvent velocity field,the equation of
motion becomes
m
i
˙
¯v
i
(t) = −mξ
i
(¯v
i
(t) − ¯u(r
i
)) +
¯
F
i
(t) +
¯
R
i
(t) (2.25)
28
2.3 Monte-Carlo procedure
A modification of the periodic boundary conditions proposed by [51] has been adopted
Figure 2.3:Sketch of shear boundary conditions for planar Couette flow.
for establishing the shear flow.The scheme used is shown in figure 2.3.In essence,the
infinite periodic system is subjected to a uniform shear in xy plane.The box above
is moving at a speed ˙γL in the positive x direction.The box below moves at a speed
˙γL in the negative x direction.When a particle in the primary box with coordinates
(x,y) crosses a horizontal border one of its images enters from the opposite side at
x

= x − ˙γtL
y
with velocity v

= v − ˙γL
y
2.3 Monte-Carlo procedure
The bonding network is itself the object of random instantaneous changes provided by
a Monte-Carlo (MC) algorithm [19] which is built according to the standard Metropolis
scheme [49].
The probability P
m,n
to go from a bounding network m to a different one n is
written as
P
m,n
= A
trial
m,n
P
acc
m,n
(2.26)
where A
trial
m,n
is the trial probability to reach a new bounding network n starting from
the old one m,within a single MC step.This trial probability is chosen here to be
symmetric as usually adopted in Metropolis Monte Carlo schemes,and is chosen to be
29
different from zero only if both bounding networks n and m differ by the status of a
single bond,say the pair of monomers (ij).To satisfy microreversibilty,the acceptance
probability P
acc
m,n
for the trial (m−→n) must be given by
P
acc
m,n
= Min[1,exp (−
(U({r},n) −U({r},m))
k
B
T
)].(2.27)
In the present case,as a single pair (ij) changes its status,the acceptance probability
takes the explicit form
P
acc
m,n
= Min[1,exp(−
(U
2
(r
ij
) −U
1
(r
ij
))
k
B
T
)] (2.28)
P
acc
m,n
= Min[1,exp(−
(U
1
(r
ij
) −U
2
(r
ij
))
k
B
T
)] (2.29)
for bond scission and bond recombination respectively.
The way to specify the trial matrix A
trial
m,n
starts by defining a range of distances,
called Γ and defined by 0.96 < r < 1.20 within the range r < R.For r ∈ Γ,a
change of bounding is allowed as long as the two restricting rules stated above are
respected.Consider the particular configuration illustrated by figure 2.4 where M=7
monomers located at the shown positions,are characterized by a connecting scheme
made explicit by representing a bounding potential by a continuous line.The dashed
line between monomers 5 and 6 represents the changing pair with distance r ∈ Γ which
is a bond U
1
(r) in configuration m but is just an ordinary intermolecular pair U
2
(r) in
configuration n.
For further purposes,we have also indicated in figure 2.4 by a dotted line all pairs
with r ∈ Γ which are potentially able to undergo a change from a non bounded state
to a bounded one in the case where the (56) pair,on which we focus,is non bounded
(state n).Note that the bond (35) is not represented by a dotted line eventhough the
distance is within the Γ range:a bond formation in that case is not allowed as it would
lead to a cyclic conformation.Note also that in state n,monomer 5 could thus form
bonds either with monomer 2 or with monomer 6 while monomer 6 can only form a
bond with monomer 5.
We now state the algorithm and come back later on the special (m−→n) transition
illustrated in figure 2.4.
During the LDdynamics,with an attempt frequency ω per armand per unit of time,
a change of the chosen arm status (bounded U
1
(r) to unbounded U
2
(r) or unbounded
30
2.3 Monte-Carlo procedure
to bounded) is tried.If it is accepted as a “trial move” of the bounding network,it
obviously implies the modification of the status of a paired arm belonging to another
monomer situated at a distance r ∈ Γ from the monomer chosen in the first place.
The trial move goes as follows:a particular arm is chosen,say arm 1 of monomer
i,and one first checks whether this arm is engaged in a bounding pair or not.
• If the selected arm is bounded to another arm (say arm 2 of monomer j) and the
distance between the two monomers lies within the interval r
ij
∈ Γ,an opening
is attempted with a probability 1/(N
i
+1),where the integer N
i
represents the
number of monomers available for bonding with monomer i,besides the monomer
j (N
i
is thus the number of monomers with at least one free arm whose distance
to the monomer carrying the originally selected arm i lies within the interval Γ,
excluding from counting the monomer j and any particular arm leading to a ring
closure).If the trial change consisting in opening the (ij) pair is refused (either
because the distance is not within the Γ range or because the opening attempt
has failed in the case N
i
> 0),the MC step is stopped without bonding network
change (This implies that the LD restarts with the (ij) pair being bonded as
before).
• If the selected arm (again arm 1 of monomer i) is free from bonding,a search is
made to detect all monomers with at least one armfree which lie in the “reactive”
distance range r ∈ Γ from the selected monomer i (Note that if monomer i is a
terminal monomer of a chain,one needs to eliminate from the list if needed,the
other terminal monomer of the same chain in order to avoid cyclic micelles config-
urations).Among the monomers of this “reactive” neighbour list,one monomer
is then selected at random with equal probability to provide an explicit trial
bonding attempt between monomer i and the particular monomer chosen from
the list.Note that if the list is empty,it means that the trial attempt to create a
new bond involving arm 1 of monomer i has failed and no change in the bonding
network will take place.
In both cases,if a trial change is proposed,it will be accepted with the probability
P
acc
m,n
defined earlier.If the change is accepted,LDwill be pursued with the newbonding
scheme (state n) while if the trial move is finally rejected,LD restarts with the original
bonding scheme corresponding to state m.
31
Figure 2.4:Exemplary configuration of a 7 monomers systemin state n where monomers
3,4 and 5 form a trimer and monomers 6 and 7 a dimer.All pairs of monomers with
mutual distances within the Γ region are indicated by a dotted or a dashed line.In the
text,we consider the Monte Carlo scheme for transitions between states n and m which
only differ by the fact that in state n and mthe 5-6 pair is respectively open or bounded.
The n →m transition corresponds to the creation of a pentamer by connecting a dimer
and a trimer while the m → n transition leads to the opposite scission.The cross
symbol on link 3-5 indicates that in state n,when looking to all monomers which could
form a new link with monomer 5,monomer 3 is excluded because it would lead to a
cyclic polymer which is not allowed within the present model.
Coming back to figure 2.4,we now show that the MC algorithm mentioned above
garantees that the matrix A
trial
mn
is symmetric,an important issue as it leads to the
micro-reversibility property when combined with the acceptance probabilities described
earlier.Let us define as P
arm
= 1/2N the probability to select a particular arm,a
uniform quantity.
If configuration m with pair (56) being “bounded” is taken as the starting con-
figuration,the number of available arms to form alternative bonds with monomer 5
and monomer 6 are respectively N
5
= 1 and N
6
= 0.Therefore,applying the MC
rules described above,the probability to get configuration n where the pair (56) has to
be unbounded is given by the sum of probabilities to arrive at this situation through
selection of the arm of monomer 5 engaged in the bond with monomer 6 or through
selection of the arm of monomer 6 engaged in a bond with monomer 5.This gives
A
trial
m,n
= P
arm

1
N
5
+1
+P
arm
1
N
6
+1
=
3
2
∗ P
arm
(2.30)
If configuration n with pair (56) being “unbounded” is taken as the starting config-
uration,the application of the MC rules lead to the probability to get configuration m
32
2.3 Monte-Carlo procedure
where the pair (56) has to be bounded is given by the sum of probabilities to arrive at
this situation through selection of the free arm of monomer 5 (which has two bounding
possibilities,namely with monomers 2 and 6) or through selection of the free arm of
monomer 6 which can only form a bond with monomer 5.This gives
A
trial
n,m
= P
arm

1
2
+P
arm
=
3
2
∗ P
arm
(2.31)
showing the required matrix symmetry.
33
34
Chapter 3
Equilibrium Properties
We have exploited the model introduced in chapter II in a series of Langevin Dynam-
ics (LD) simulations at equilibrium with different bonding energy parameter W and
different number density φ,and hence,in this chapter,we will present the resulting
equilibrium static and dynamic properties of cylindrical micelles.
Within the static properties,the theoretical prediction of the chain size distribution
has been given by Cates [2] and tested in great detail by Monte Carlo simulation by
Wittmer et al.[13,52] who also investigated conformational properties of chains at
equilibrium.For static properties,our aim is thus mainly to check that our results on
a different model are compatible with the analysis of the previous works [13,52].
The two main parameters governing the static properties are the monomer density
φ and the end cap energy E.Our choice of these two parameters is set up with the
aim to simulate three thermodynamic states corresponding to a dilute solution and two
semi-dilute solutions at the same φ but different E,leading to the system with two
average chain lengths ∼ 56 and ∼ 150.
For static properties,the distribution of chain lengths,the gyration radius,and
the end-to-end distance versus chain length will be analyzed and compared with pre-
vious studies.We have also investigated g
ee
(r) the pair correlation function of end
monomers and P(r) the distribution of bond lengths which are quantitative pertinent
in the microscopic formulation of the kinetic rate constants.
The dynamic properties of the micelles which are given and interpreted in this
chapter form the core of the work.The detailed trajectory of diffusing micelles which
are continuously breaking and recombing allows us to analyze the microscopic origin of
rate constants in terms of structural features (e.g.chain end pair correlation function)
35
and dynamic quantities to be related to the statistics of life times of a newly created
chain end.For the latter,we show that the Poisson statistics dominates at long times
while a fraction of correlated recombinations happen at short times.Exploiting these
microscopic features,we characterize the macroscopic scission energy E and the barrier
of recombination B and estimate their values for various state points investigated.
Some macroscopic dynamic properties are then studied with an accent on the mod-
ification of various dynamic relaxation processes due to the scission-recombination pro-
cess.We investigate the monomer diffusion and the stress relaxation function.Finally,
we perform a T-jump experiment in order to point out that the previously estimated
macroscopic kinetic constants are indeed the key parameters governing the relaxation
of the chain length distribution.
3.1 Static properties
The main aim of this section is to test the model of chapter II by comparing the struc-
tural properties,including the average chain length,the distribution of chain lengths
and the conformations of the chains,with the theory [9,27] and previous simulation
works [13,52].
3.1.1 List of simulation experiments and chain length distribution
The model is studied at three state points.The number of monomers,the number
density φ,the energy parameter W,and the attempt frequency ω per arm are chosen
for
1.A solution at the number density φ = 0.05 and an energy parameter W = 8.
The number of monomers is M=1000.The attempt frequencies of bond scis-
sion/recombination per arm ω are 0.1,0.5,1 and 5.This choice will be shown to
lead to a dilute solution.
2.A solution at the number density φ = 0.15 and an energy parameter W = 10.
The number of monomers is M=1000.The attempt frequencies of bond scis-
sion/recombination per arm ω are 0.1,0.5,1 and 5.(will be shown to be a
semi-dilute solution)
36
3.1 Static properties
3.A solution at the number density φ = 0.15 and an energy parameter W = 12.
The number of monomers is M=5000.The attempt frequencies of bond scis-
sion/recombination per arm ω are 0.02,0.06,0.1,0.5 and 1 (will be shown to be
a semi-dilute solution).
Each system evolves according to the Langevin Dynamics algorithm with time step
Δt = 0.005 and is subject to randomtrials of bond scission/recombination with the arm
attempt frequency ω.All the experiments and the results of static properties which
include the average chain length L
0
,the mean square end-to-end vector,the radius
of gyration and L
0
/L

,the number of blobs in a chain of length L
0
,are listed in the
table 3.1.Where


R
2

and


R
2
g

are defined as


R
2

L
0
=
L
0
−1
X
n=1
L
0
−1
X
m=1
h¯r
n
 ¯r
m
i (3.1)


R
2
g

L
0
=
1
2L
2
0
L
0
X
n=1
L
0
X
m=1


(
¯
R
n

¯
R
m
)
2

,(3.2)
where ¯r
n
=
¯
R
n+1

¯
R
n
.And L

is defined by equation (1.21) (also see section 3.1.2).
As shown in the table,we observed that all the static properties are independent
of ω and therefore,all data can be averaged over all ω values.
3.1.1.1 The dilute case
From Table 3.1,it can be observed that the first state point experiment (W = 8,φ =
0.05) is a dilute solution since its average chain length L
0
= 11.48(1) is much smaller
than its crossover value at the monomer number density as calculated by equation (1.20),
L

= 50.5.Dilute solution conditions are confirmed by the chain length distribu-
tion shown in figure 3.1.For dilute conditions,a distribution given by (1.30) is ex-
pected [13,52].We fit our data with a single parameter B
1
of function (1.30),where
L
0
is given its computed average value and where γ is given its expected value,γ = 1.165
[13].The fit gives B
1
= 1.08.This curve is significantly better than the simple expo-
nential distribution expected for ideal or semi-dilute chain.If γ is left as a second free
parameter in the fit,it gives γ = 1.161 which is also very close to its expected value.
37
Table 3.1:List of simulation experiments and the values of static properties.W is the
scission energy parameter,φ the number density,T
s
the total simulation time,ω the
attempt frequency,L
0
the average chain length,
q
hR
2
i
L
0
the end-to-end distance of
the average chain,and
q


R
2
g

L
0
its radius gyration.L
0
/L

is the ratio of the average
chain length over the blob length L

W
φ
T
s
ω
L
0
q
hR
2
i
L
0
q


R
2
g

L
0
L
0
/L

8
0.05
2.5 ∗ 10
5
0.1
11.48(6)
4.93(2)
1.9(2)
0.2
8
0.05
2.5 ∗ 10
5
0.5
11.46(2)
4.92(1)
1.9(1)
0.2
8
0.05
2.5 ∗ 10
5
1.
11.50(3)
4.94(1)
1.9(1)
0.2
8
0.05
2.5 ∗ 10
5
5.
11.51(2)
4.95(1)
1.9(1)
0.2
10
0.15
2.5 ∗ 10
5
0.1
56.2(6)
12.2(1)
4.8(2)
4.7
10
0.15
2.5 ∗ 10
5
0.5
56.2(1)
12.4(1)
4.9(3)
4.7
10
0.15
2.5 ∗ 10
5
1.
56.2(2)
12.3(1)
4.9(3)
4.7
10
0.15
2.5 ∗ 10
5
5.
56.6(2)
12.4(1)
4.9(4)
4.7
12
0.15
6.25 ∗ 10
5
0.02
151(4)
20.9(5)
8.3(4)
12.6
12
0.15
4.5 ∗ 10
5
0.06
153(4)
21.5(7)
8.5(5)
12.6
12
0.15
4.5 ∗ 10
5
0.1
150.7(5)
21.2(2)
8.4(3)
12.6
12
0.15
3 ∗ 10
5
0.5
150.3(5)
20.9(1)
8.4(2)
12.5
12
0.15
3 ∗ 10
5
1
150.4(6)
20.92(3)
8.5(4)
12.5
3.1.1.2 The semi-dilute case
Both the second state point (W = 10,φ = 0.15) and the third state point (W =
12,φ = 0.15) experiments are found to be in semi-dilute regime,since their average
chain lengths are 56.4(1) and 151.4(4),respectively,which are several times larger
than their crossover value L

≈ 12 at φ = 0.15 according to equation (1.20).Semi-
dilute solution conditions are confirmed by the observation of a simple exponential
distribution of chain lengths.In figure 3.2 and figure 3.3,we show,for the second
and the third state point experiments respectively,the agreement of our data with the
expected simple exponential distribution (1.34).
3.1.2 Chain length conformational analysis
In this subsection we are interested in the conformational properties of micelles in
equilibrium and studied as a function of chain size within our polydisperse system.
We have calculated the mean square end-to-end distance and the mean square radius
of gyration < R
2
> (L) and < R
2
g
> (L) averaged over subsets of chains of length L.
38
3.1 Static properties
0 1 2 3 4 5 6 7 8 9 10
L/L
0
10
-8
10
-7
10
-6
10
-5
10
-4
10
-3
c0(L)
Figure 3.1:Distribution of chain lengths for the dilute case in good solvent.The data
(squares) are fitted using equation (1.30) with the single parameter B
1
and imposed
value γ = 1.165 (continuous curve).The fit gives B = 1.08.The dashed line shows
the simple exponential distribution c
0
(L) ∝ exp(−L/L
0
) which does not fit the data.
Figure 3.4 and figure 3.5 are the results for the dilute state point and the two semi-dilute
state points experiments respectively.As shown in figure 3.4 relative to the dilute case,
the squares and circles represent our data of < R
2
> and < R
2
g
> respectively.With
the solid line and the dashed line,we indicate,for long chains (L ￿ 25),standard
power law scaling L

with ν = 0.588.The data show a reasonable agreement of this
asymptotic regime in the range of application.
The chains in the semi-dilute system are expected to behave as ideal chains for
L ≫ L

≈ 12 where L

is estimated from equation.(1.20) with φ = 0.15.Figure 3.5
shows the


R
2

and


R
2
g

of the two semi-dilute cases.Results of


R
2

and


R
2
g

for the two cases,are superimposed.We indicate for the long chains the ideal chain
conformation


R
2

and


R
2
g

∝ L.The agreement of the fitting lines with simulation
data appears to start at L ￿ 60.
The conformational properties of the chain in the two semi-dilute cases are found
39
0 1 2 3 4 5 6 7
L/L
0
10
-8
10
-7
10
-6
10
-5
10
-4
c0(L)
Figure 3.2:Distribution of chain lengths for the state point W = 10,φ = 0.15 (ac-
cumulated over all ω values).The data (circles) are fitted very well by the simple
exponential function,equation (1.34),with imposed average chain length L
0
= 56.4.
The Dashed line shows the function exp(−1.165L/L
0
) which does not fit the data.
to be identical as expected.The size R of a chain of length L > L

is
R = bL

ν

L
L


0.5
(3.3)
where b is the monomer size and ν is the good solvent scaling exponent while L

(function of φ only) is the same in both cases.
3.1.3 Pair correlation function of chain ends and the distribution of
bond length
Equilibrium polymers are polydisperse polymers endowed with scission and recombi-
nation processes.Whereas a scission can happen between any bounded pair,a recom-
bination may happen only with two chain ends for linear chains.Thus it is interest-
ing to study the spatial distribution of the chain ends and the bond length distribu-
tion.For the dilute case the population of bonds ready to open within the Γ range
(−0.96 < r < 1.2) is < N

>= 552 and represents 61% of the bonded pairs < N
1
>,
while the population of free arm pairs ready to close,again within the same Γ range,is
40
3.1 Static properties
0 1 2 3 4 5 6
L/L
0
10
-8
10
-7
10
-6
10
-5
c0(L)
Figure 3.3:Distribution of chain lengths for the state point W = 12,φ = 0.15 (ac-
cumulated over all ω values).The data (circles) are fitted very well by the simple
exponential function,equation 1.34,with imposed average chain length L
0
= 151.4.
The Dashed line shows the function exp(−1.165L/L