Computer Simulation Studies of Chain Dynamics in Polymer

Thus, for ϕ ReRe (i) the concept of the initial relaxation time is not very useful in practice. However, with respect to the correlator describing th...
0 downloads 0 Views 4MB Size
Article pubs.acs.org/Macromolecules

Computer Simulation Studies of Chain Dynamics in Polymer Brushes Daniel Reith,† Andrey Milchev,‡ Peter Virnau,† and Kurt Binder*,† †

Institut für Physik, Johannes Gutenberg-Universität Mainz, Staudinger Weg 7, D-55099 Mainz, Germany Institute of Physical Chemistry, Bulgarian Academy of Sciences, 1113 Sofia, Bulgaria



ABSTRACT: Center-of-mass and single monomer motion in grafted chains comprising a strongly stretched polymer brush in thermal equilibrium are studied by large scale molecular dynamics and Monte Carlo simulations of a coarsegrained model. Good solvent conditions are assumed. Our findings seriously question earlier theoretical predictions about the relaxation described by Rouse dynamics of brush coatings. Thus, the correlation functions of parallel and perpendicular components of the mean distance of the center-of-mass from the grafting site, the squared gyration radius and end-to-end distance, are found to deviate strongly from a simple exponential decay. While the relaxation times extracted from the initial slope of the center-of-mass and the gyration radius correlation functions are roughly compatible with the theoretical prediction τ(i) z ∝ 2 N3, τ(i) xy ∝ N , the relaxation times extracted from the late exponential decay scale as τ(z f) ∝ τ(xyf) ∝ NΔ with an effective exponent Δ ≈ 3.7 for polymers of length 64 ≤ N ≤ 256. Moreover, by studying the dynamics of individual monomers from the grafting site (i = 1) up to the free chain end (i = N), it is found that monomers near the middle (i = N/2 + 1) exhibit the slowest relaxation, although the range of their mean square displacements is clearly much smaller than that of the end monomers. Similar results are found also in the case of brushes formed from ring polymers.

τz ∝ ⟨(δz)2 ⟩/Ddumb ∝ a 2(σga 2)2/3 N3ζ0

I. INTRODUCTION Recently, there has been much work devoted to the understanding of polymer brushes, dense polymeric layers of flexible polymers grafted by a special chemical group to a flat or curved substrate.1−9 Both the static properties and the response to shear and other out of equilibrium conditions have been intensively studied, see9−13 for recent examples and further references. However, much less attention has been paid to the monomeric and chain motions in polymer brushes in equilibrium.14−18 Of course, understanding the dynamics of fluctuations of a system near equilibrium normally is a prerequisite to fully understand the far from equilibrium behavior.5,9,11−13 Thus, the present paper attempts some steps to close this gap. We will not include solvent molecules explicitly, and thus the dynamics of comparably short chains in solution under the considered good solvent conditions corresponds to the Rouse- rather than the Zimm model, lacking a proper description of hydrodynamic interactions in the dilute limit.19 The polymer brushes to be considered here, however, imply rather large grafting densities, and hence the conditions inside the brush correspond to at least semidilute or even concentrated solutions, where hydrodynamic interactions are screened out to a large extent.19 Thus, we feel that our model should capture the actual behavior at a qualitative level. The existing theory14 is based on an even simpler concept, where the relaxation of the polymer chain is described by a simple dumbbell model.20 For free-draining chains with N monomeric units the friction coefficient then is ζ = ζ0N, ζ0 being the segmental friction coefficient. The diffusion coefficient then is Ddumb = kBT/ζ = kBTζ0−1N−1, and the relaxation time of the z component is14 © 2012 American Chemical Society

(1)

a being the size of a monomer and σg the grafting density. Here we have used the fact, shown by the self-consistent field theory (SCFT) of polymer brushes21−23 and confirmed by simulations,24−29 that the mean-square fluctuations of the end-toend distance of the dumbbell ⟨(δz)2⟩ scale with chain length N and grafting density σg in the same way as the square of the brush height h, ⟨(δz)2⟩ ∝ h2, with h ∝ a(σga 2)1/3 N

(2)

In eqs 1, 2 a is the linear dimension of a monomeric unit (so that at the Theta point we would have for the end-to-end distance30 Re ∝ a√N). Evidently due to this anomalous large fluctuation14 ⟨(δz)2⟩ ∝ 2 N eq 1 differs from the standard scaling of the Rouse model (with Flory exponent ν = 1/2)19 τRouse ∝ Re 2/Ddumb ∝ a 2N 2ζ0

(3) 2

However, a scaling of the relaxation time with N , as in eq 3, is predicted if we consider the fluctuations of the end-to-end vector of the dumbbell in the x- or y directions parallel to the grafting surface, for which one predicts 15,26 ⟨(δx)2 ⟩ = ⟨(δy)2 ⟩ ∝ a 2(σga 2)−1/6 N

(4)

Received: December 20, 2011 Revised: March 1, 2012 Published: May 10, 2012 4381

dx.doi.org/10.1021/ma202745b | Macromolecules 2012, 45, 4381−4393

Macromolecules

Article

with constants k, R0 chosen as k = 30, R0 = 1.5. The grafting of the first (mobile) monomer is effected by the same potential, eq 8. The model used for the MC simulations was chosen somewhat different, for two reasons: (i) the model has been extensively studied in previous work both in bulk melts and in solutions, and hence is well characterized34−36 (ii) Because of its particular short-range of intermonomer interaction, a rather fast simulation code could be written. An additional bonus is that by comparing results from two somewhat different coarsegrained models we can better tell which properties are still model-dependent and which properties describe universal features. Of course, only the latter are of particular interest here. In this model, bonded monomers interact with a different type of FENE potential, namely35

and using the same type of scaling argument for the relaxation time as used in eq 1, one finds 15 τx ∝ ⟨(δx)2 ⟩/Ddumb ∝ a 2(σga 2)−1/6 N 2ζ0

(5)

15

Early Monte Carlo studies reported estimates compatible with eqs 1, 5, but the present work will show that the situation is not so simple: we will show that τx and τz actually have the same N-dependence, and the effective exponent Δ in the relation τx ∝ τz ∝ NΔ actually exceeds the prediction Δ = 3 (eq 1) significantly. This finding is surprising, since eq 1 is also corroborated by a more elaborate treatment, extending the Rouse model to a chain in a polymer brush where the neighboring chains are described in terms of a self-consistent potential.17 The outline of our paper is as follows. In section II, we briefly describe the model and simulation methodology. In section III, we describe the pertinent static results, to show that our chain lengths are large enough to verify the scaling relations for the static linear dimensions (eqs 2, 4). We shall also include the case of brushes formed from ring polymers,31 for which a relation ⟨(Rgx)2⟩ = ⟨(Rgy)2⟩ ∝ N0.8 instead of eq 4 was established.31 Section IV describes our findings on the relaxation times, and section V summarizes some conclusions.

U ′FENE (S) = −U0 ln[1 − (S − S0)2 /(S0 − Smin)2 ]

where the bond length (S ) can vary between S min < S < S max, where S max = 1 is chosen as the unit of length, S 0 = 0.7, S min = 0.4. The constant U0 is chosen as U0 = 20. Nonbonded interactions are described by a Morse potential34−36 UM(r ) = εM (exp[−2α(r − rmin)] − 2 exp[−α(r − rmin)]) (10)

with parameters εM = 1, rmin = 0.8, α = 24. Because of the attractive character of this potential at small distances, under dilute solution conditions for N → ∞ a well-defined Θ point occurs, which has been estimated to be35 kBTθ/εM ≈ 0.62. However, here we study the model only for kBT/εM = 1, so we consider only the good solvent regime, just as for the athermal model defined by eqs 7, 8 that is used in the MD work. A possible extension of the study to the regime of Θ solvents or poor solvent conditions is left to future work. Dynamics is associated with Monte Carlo simulations as usual,36,37 via the Metropolis algorithm: a time step consists in the random selection of a monomeric unit in the system for which one attempts a random move from its old position ri⃗ to a new position r′⃗ i = ri⃗ + Δr.⃗ The displacement Δr ⃗ is randomly chosen from a cube of linear dimension δ centered at the old position, with δ = 0.5. The energy change in the system is calculated {using eqs 9 and 10} and the attempted move is actually carried out (or not), subject to the standard Metropolis test.36,37 Note that in the present case the displacement scale δ = 0.5, results in a typical acceptance rate 0.14 of attempted moves. The resulting reduced time unit of the MC simulation then is 1 MCS when every monomer in the system on average has had the opportunity to move once. As is well-known, this algorithm for chains in dilute solution leads to Rouse-like relaxation of the polymer chains;34−37 hydrodynamic effects mediated by solvent molecules, that are not explicitly considered, are missing, of course. Since the potentials eqs 9 and 10 are chosen such, that in the course of these random motions of the monomeric units a move where two chains cross through each other cannot occur, the model fully respects entanglement constraints on the chain dynamics.36,37 The molecular dynamics simulations are based on the use of the Velocity Verlet algorithm for the integration of Newton’s equations of motion for the particles.2,5,33 The mass m of the monomeric unit is chosen as m = 1, so the MD time unit is τ = σ(ε/m)1/2 = 1 as well. MD runs were carried out on graphics processing units (GPUs) using the HooMD-blue r4001 code.38 Note that due to the use of GTX 480 CPUs, a speed-up of roughly 70 in comparison to a run on a single i7 CPU core for

II. MODELS AND SIMULATION METHODS The substrate surface for our polymer brushes is taken to be a piece of the xy-plane with linear dimensions L × L and periodic boundary conditions. In the molecular dynamics (MD) simulations, the wall is represented by particles fixed on the sites of a square lattice with lattice constant a0 and interacting with the effective monomers with a potential of the Weeks− Chandler−Andersen (WCA)32 type ⎛ 4ε [(σ /r )12 − (σ /r )6 + 1/4], r ≤ 21/6σ w w w w VWCA(r ) = ⎜ ⎜ 1/6 r > 2 σw ⎝ 0, (6)

where εw and σw characterize the strength and range of this repulsion. For the Monte Carlo (MC) simulation, instead a simple hard wall {VHW(z ≤ 0) = ∞, VHW(z > 0) = 0} is assumed. These potentials are impenetrable for the monomeric units, whose z-coordinates are hence restricted to the positive half space, z > 0. For most of our MD simulations chains are grafted at a density σg = 0.125σ−2, σ being the range parameter of the WCA potential that acts between effective monomers ⎡ 1⎤ UWCA(r ) = 4ε⎢(σ /r )12 − (σ /r )6 + ⎥ , ⎣ 4⎦

r ≤ rc = 21/6σ (7)

while UWCA(r ≥ rc) = 0. For simplicity, we take εw = ε = 1 as the unit of energy (choosing also the temperature T = 1.0ε/kB = 1, with Boltzmann’s constant kB1), and σw = σ = 1 is our unit of length for the MD simulations. The lattice constant a0 of the substrate in the MD model was then chosen as a0 = σ/2 = 1/2 in order to avoid interpenetration of the substrate by the brush monomers. As is standard,2,5,33 bonded monomers are connected by additional anharmonic “springs” described by the finitely extensible nonlinear elastic (FENE) potential UFENE(r ) = −0.5kR 0 2 ln[1 − (r /R 0)2 ],

r < R0

(9)

(8) 4382

dx.doi.org/10.1021/ma202745b | Macromolecules 2012, 45, 4381−4393

Macromolecules

Article

Figure 1. Log−log plot of (a) ⟨Rez⟩, ⟨Rgz2⟩1/2 and (b) ⟨Rez2⟩ − ⟨Rez⟩2, ⟨Rgxy2⟩ versus chain length NL for open chains and NR/2 for rings, as resulting from the MD simulation. Straight lines indicate power law fits, slopes are shown to indicate effective exponents. Note that in part a data for open chains and corresponding data for rings (NR/2 = N) coincide within the statistical errors, which do not exceed the size of the symbols, so only two rather than four straight lines are shown. Inset in part b shows a plot of the reduced ratio ⟨Rez2⟩/⟨Rez⟩2 − 1, to illustrate the slow approach to the thermodynamic limit.

never occur, the fact that no permanently entangled pair of rings can occur is strictly preserved during our simulation. This system hence is fundamentally different from the so-called “loop-brushes”,41 where linear chains are grafted with both ends at the substrate at randomly chosen sites, leading to a brush of loops with many permanent entanglements, and hence the substrate is covered with a kind of soft polymer network. The reason for considering a brush formed from strictly noncatenated rings is that the linear dimensions of these ring polymers have anomalous properties.31 Similar to melts of noncatenated ring polymers, where the gyration radii Rg are known to scale as42−51 Rg ∝ N νeff with an (effective) exponent 1/3 ≤ νeff ≤ 0.4 we have found that the parallel component of the gyration radius, (⟨Rgxy2⟩)1/2, shows such a scaling in ring polymer brushes31

the present application is obtained. This substantial gain in speed was absolutely crucial for being able to sample autocorrelation functions of chain linear dimensions down to small values (of order 10−2) with the required accuracy. A critical issue for any MD simulation is the proper thermalization. To achieve a fast approach to equilibrium, and sample equilibrium properties accurately, the use of standard dissipative particle dynamics (DPD) thermostat39 with parameters γ = 0.5 (friction coefficient), cutoff for the interaction range rcut = 1.25 × 21/6 and integration time step δt = 0.005 was found useful. However, for recording autocorrelation functions in time this thermostat is less useful, since collective density oscillations in the brush occur and are only rather weakly damped (in a real system, the interactions of the solvent particles with the effective monomers, together with the suppression of hydrodynamic effects close to the grafting surface by interactions of solvent particles with the walls, are expected to remove these effects to a large extent). Therefore, we have found it more useful to use again the standard Langevin thermostat2,5,33 with friction coefficient γ = 0.5 which strongly damps out any oscillatory motions (but is known to remove also hydrodynamic interaction in MD simulations of polymer solutions40). For the MD runs, L was adjusted such that we have either 5 = 102400 or 5 = 56200 effective monomers in the system, and chain lengths N = 16,32,64,128,256 and 512 were studied. Thus, the linear dimensions in the x, y directions are very large and experience with previous works on related models (where distinctly smaller linear dimensions were used) ensures that finite size effects need not be considered. For the longest chains, an equilibration time of τeq ≈ 106 MD time units was needed. For the MC simulations, somewhat smaller systems were studied, containing 5 = 32768 effective monomers in total. Equilibration for the longest chains needed τeq ≈ 106MCS, and averages were taken carrying out MC runs over typically 6 ÷ 10 × 106 MCS. In addition to linear chains of length N, also ring polymers of length NR = 2N at grafting density σg/2 were studied, such that the total number of monomers was strictly the same in both cases. The preparation of the initial state was always done such that no catenation (permanent entanglement) between any pair of rings whatsoever did occur. Since both the MD algorithm (due to the choice of the potentials {eq 7, 8}) and the MC algorithm are constructed such that in the course of random motion of the monomers any chain intersection can practically

ν

⟨R gxy 2⟩1/2 ∝ NReff ,

νeff ≈ 0.40,

32 ≤ NR ≤ 512

(11)

In the spirit of eq 5, one hence would predict that τxy ∝ ⟨Rgxy2⟩/Ddumb ∝ N1−2νeff ≈ N1.8. Checking this prediction will provide one further test of the concepts outlined in the Introduction (eqs 1−5).

III. SIMULATION RESULTS FOR STATIC BRUSH PROPERTIES As described in the Introduction, the current understanding of chain dynamics in polymer brushes rests on a somewhat qualitative extension of the Rouse model to the case of such stretched chains. It is assumed that the chains are long enough that simple power laws hold21−23 ⟨Rez⟩ ∝ N , ⟨Rez 2⟩ − ⟨Rez⟩2 ∝ N 2

(12)

⟨R gz 2⟩ ∝ N 2

(13)

⟨Rexy 2⟩ ∝ ⟨R gxy 2⟩ ∝ N

(14)

Note that in the present work we are not concerned with a systematic study of the dependence of relaxation times on the grafting density σg. Equations 1 and 5 consider only the dependence of the relaxation time on grafting density due to the fact that the linear dimensions depend on grafting density. A possible dependence of the friction coefficient ζ0 on grafting density is not considered. However, in view of the fact that in the analogous problem of Rouse dynamics and its dependence on density in concentrated solutions where ζ0 is found to 4383

dx.doi.org/10.1021/ma202745b | Macromolecules 2012, 45, 4381−4393

Macromolecules

Article

Figure 2. (a) Probability distribution Pi(x) of monomer i (label along the chain contour) for a chain of length N = 64, for several choices of i, namely i = 16, 32, 48 and the chain end i = N. (b) Same as part a, but for the z-coordinate of monomer i, Pi(z) versus z. All data are from the MD simulation.

depend strongly on concentration,52,53 a study of the dependence of ζ0 on σg is warranted, but is left for future work. For a study of this problem in a related model with only two chain lengths N = 64 and N = 128, see ref 18. Figure 1 presents our data on these chain dimensions in polymer brushes on a log−log plot. Estimates for rings as well as linear chains are included in this and some of the following figures. We will refer to chain length N as N = NL in case of linear chains and N = NR/2 in case of rings to avoid any confusion. While the first relation of eq 12, as well as eqs 13 and 14, is verified nicely, one can see that in the range of chain lengths N that is studied, the second relation of eq 12 is not yet fully verified over the whole range of chain lengths N. There rather seems to be a somewhat slow approach to the scaling limit, cf. inset of Figure 1b. This slow approach is also obvious by the fact that there is clear curvature visible on the log− log-plot in Figure 1b when smaller values of N are included. For the rings, however, we observe ⟨Rex2 ⟩ = ⟨Rey 2⟩ ∝ NR 0.95 ,

32 ≤ NR /2 ≤ 512

grafted ring corresponds to two grafted linear chains, and then the total monomer distribution ρ(z) in such brushes is found to be indistinguishable in both cases31). Figure 3 shows the probability distribution of the z component Rgz of the radius of gyration.

(15)

Figure 3. Probability distribution P(Rgz) of the z component Rgz of the gyration radius of polymers, comparing rings with NR monomers with linear chains with N = NL = NR/2 monomers, but twice the grafting density. NR = 128, 256, 512, and 1028; as indicated (lower part). Upper part shows corresponding scaling plots, P(Rgz/⟨Rgz⟩) versus Rgz/⟨Rgz⟩, for rings (left) and open chains (right).

for the transverse components of the end-to-end distance. The fact that the fluctuation of the end monomer distance from the grafting surface shows a slow approach to the scaling limit is plausible when one studies the distribution functions of the positions of individual monomers in a brush (Figure 2). While the distribution of the x- or y-coordinates are simple gaussians centered around zero, and the widths of this Gaussian increases steadily as one progresses from the monomers near the grafting site to the free end (Figure 2a), the distribution of the z-coordinate has a complicated asymmetric shape (Figure 2b), but still differs substantially from the asymptotic form derived by the self-consistent field theory,21−23 which predicts a square-root singularity near the brush height, P(z) ∝ (h − z)1/2. Nevertheless, Figure 2b) qualitatively confirms the prediction21−23 that the free chain end is not localized in the “last blob” (as anticipated by the Alexander−de Gennes model54,55 of brushes), but can occur anywhere throughout the brush. The same conclusion did emerge from the end-monomer distribution function ρe(z) recorded for the corresponding model by MC simulation, that was already presented in ref 31, and hence shall not be reproduced here. We now turn to a comparison of static chain properties of ring polymer brushes with corresponding linear polymer brushes (recall that NR = 2N and σg,r = σg/2, so that every

One sees that with respect to the z component of the gyration radius, ring polymer brushes and the corresponding brushes formed from linear chains again are indistinguishable. But the deviations from a good “data collapse” near the peaks in the corresponding scaling plots (upper part of Figure 3) indicate again a rather slow approach to the asymptotic scaling limit, as pointed out already for the end-to-end distance. The story is different, however, for the xy components of the gyration radius (Figure 4). We see very good “data collapse” in the scaled distributions P(Rgx) vs Rgx/⟨Rgx⟩, both for rings (left upper part) and for open chains (right upper part), for the considered range of N or NR, respectively. But the shape of the ring distribution differs clearly from the shape of the chain distribution, as highlighted by the lower part of the figure. For the chains the distributions clearly are wider. Of course, within statistical errors (which do not exceed the size of the symbols in Figure 4) distributions P(Rgx) and P(Rgy) coincide since x and y directions are equivalent. 4384

dx.doi.org/10.1021/ma202745b | Macromolecules 2012, 45, 4381−4393

Macromolecules

Article

⟨Rex2⟩, ⟨Rey2⟩ for linear chains and rings disagree. This difference was expected from the corresponding difference in gyration radii. However, while we did find ⟨Rgxy2⟩ ∝ N0.8 for the rings, we now find ⟨Rex2⟩ ∝ N0.95 (cf. eq 15). This difference in (effective) exponents between ⟨Rexy2⟩ and ⟨Rgxy2⟩ for rings clearly is unexpected, and it is of course doubtful whether or not the asymptotic region was really reached. But since the structure of a polymer in a brush is rather heterogeneous, displacements of monomers near the outer region of the brush can exceed those in the bulk of the brush, there is not necessarily a contradiction between eqs 11 and 15. For the linear chains, no such pattern occurs, ⟨Rex2⟩ ∝ N as well as ⟨Rgxy2⟩ ∝ N, as expected. Figure 4. Probability distribution P(Rgxy) of the xy component Rgxy of the gyration radius of polymers, comparing rings with NR monomers with linear chains with N = NL = NR/2 monomers, but twice the grafting density. Several values of NR are shown, as indicated. The upper part shows the corresponding scaling plots.

IV. RELAXATION TIMES A standard tool to investigate relaxation of some dynamical variable A(t) is to consider its normalized autocorrelation function ϕAA(t ) = [⟨A(0)A(t )⟩ − ⟨A2 ⟩]/[⟨A2 ⟩ − ⟨A⟩2 ]

Figure 5 finally shows log−log plots of the average mean square gyration radii in the z direction and in the xy direction, comparing linear chains and rings, as obtained from both MD and MC runs. As we have discussed already in our preliminary communication,31 the component parallel to the grafting plane satisfies eq 11 while the data for ⟨Rgz2⟩ of rings coincide precisely with those of the linear chains. In order to discuss a quantity for ring polymers that is analogous to the end-to-end distance Rek(k = x,y,z) of linear polymers, we can take the coordinate of the monomer with index i = NR/2 + 1 (i = 1 being the monomer of the ring which is connected to the substrate via a spring, and monomers are consecutively labeled along the contour, i = 1, 2, 3, ..., NR being again a nearest neighbor of the first monomer). Figure 1a shows a beautiful linear scaling of ⟨Rez⟩ vs N defined for rings in this way, and also for this quantity data for rings and linear chains superimpose precisely. For ⟨Rez2⟩ − ⟨Rez⟩2, surprisingly, the data superimpose only for very large N, for smaller N the data for rings fall systematically below those of the corresponding linear chains: but this has the effect that the second eq 12 works very well for the rings, ⟨Rez2⟩ − ⟨Rez⟩2 ∝ N2 describes all data for N ≥ 32 nicely, unlike the corresponding data for linear chains (Figure 1b). The transverse components

(16)

Note that in thermal equilibrium, any time-displaced correlation function ⟨A(t1)A(t2)⟩ (with t2 > t1) must show time translation invariance, ⟨A(t1)A(t2)⟩ = ⟨A(0)A(t2 − t1)⟩, and for t2 − t1 → ∞ correlations must decay, so ⟨A(0)A(t2 − t1)⟩ → ⟨A(0)⟩⟨A(t)⟩ = ⟨A⟩2, observables ⟨A(t)⟩ being independent of time. Then often it is true that the asymptotic decay of ϕAA(t) toward zero is exponential, ϕAA(t ) ∝ exp( −t /τAA), t → ∞

(17)

and such an assumption in fact is compatible with the relaxation of various polymer brush autocorrelation functions (Figure 6a). Since ϕAA(t) as defined in eq 16 starts out at ϕAA(t = 0) = 1, one often assumes that eq 17 essentially holds at all times, and (e) then a time τ(e) AA is defined from the condition that ϕAA(t = τAA) has decayed to 1/e. Obviously, such an approach (although it has been used in the past15) is clearly inappropriate here. Since ϕAA(t) is initially strongly curved, the regime where eq 17 holds is only reached after ϕAA(t) has decayed to values smaller than 10−1. In an early work on this problem,15 only rather limited “Monte Carlo data” were available in the regime ϕAA(t) ≥ 0.2, and hence this problem was not noticed. The fact that ϕ AA(t) does not just decay with a single relaxation time τ AA,

Figure 5. Mean squared components of the gyration radii of ring polymers and equivalent linear polymers of chain length N (N = NR/2, σg = 2σg,r) shown as log−log plot versus N including both xy components (Rgxy) and the z component (Rgz2) Part a refers to the MD model and part b to the MC model. 4385

dx.doi.org/10.1021/ma202745b | Macromolecules 2012, 45, 4381−4393

Macromolecules

Article

Figure 6. (a) Semilog plot of the autocorrelation functions ϕAA(t) versus time for the choices A = R2ez, R2ex, R2gz, and R2gxy, as obtained from MD. (b) Autocorrelation functions parallel (||) and perpendicular (⊥) to the grafting plane for the gyration radius R2g −, and the center-of-mass CM−motion in a CM polymer brush with N = 64 and σg = 0.125 from MC simulation data. Inset shows scaling of τCM xy , τz , with N. Part c shows an effective time-dependent (MD) relaxation time τefAAf (t), where the exponential decay of 50 neighboring data points (from t − 25 to t + 25) was fitted in terms of eq 18. The effective time-dependent (MD) relaxation time for t < 200 in part d was fitted in terms of eq 18, too, but only for two neighboring points (from t to t + 1). Dashed horizontal lines denote the estimated initial slope of the relaxation. Note that the data of parts c and d refer to the center-of-mass motion, for which the Rouse model predicts the same N-dependence for the initial and the asymptotic relaxation time (see Appendix).

in the coordinate system of the corresponding grafting site. Figure 6c,d show that τAA may exceed τ(i) AA by about an order of magnitude. Only when the curves in Figure 6c have settled down to horizontal plateaus can one be sure that the asymptotic region has been reached (which in practice is difficult to judge, since the statistical fluctuations at late times, may prevent a clear-cut conclusion, see Figure 6c). A small value τ(i) AA at the initial relaxation may still rise to much larger values at late times. Figure 6 just presents typical examples for a relatively short chain (N = 64); data for larger N presumably are even more limited, since the largest relaxation time increases at least like N3.7 with N, see below. Figure 7a presents our results for all the relaxation times that were estimated from our simulations (due to lack of statistics in the MC simulation at very large times, only initial relaxation times extracted from MC are shown, while Figure 7a provides a tentative estimate of asymptotic relaxation times estimated via MD, as shown in Figure 6). It is seen that our MD data for N ≥ 64 are compatible with a power law with an empirical exponent Δ,

describing the asymptotic decay, but rather with a full spectrum of relaxation times, is well-known in the context of the Rouse model,19 of course. It also becomes evident when we define an effective time-dependent relaxation time τeff AA(t) as follows ⎛ d[ln ϕ (t )] ⎞−1 eff AA τAA (t ) = − ⎜ ⎟ d t ⎝ ⎠

(18)

Of course, eq 18 is most useful when τefAAf (t) is of the same order of magnitude for all times: in a polymer context we should hence consider quantities for which it is expected that ef f even the initial relaxation time τ(i) AA = limt→0τAA(t) scales with the same power of N as the asymptotic relaxation time τAA, defined in eq 17. It is important to note that this is not the case if A is the end-to-end vector R⃗ e, but it is expected to be true when A is related to the center-of-mass motion (for a more detailed discussion of the problem see Appendix). The simplest autocorrelation functions to consider are those relating to the center-of-mass motion, ϕRxCM RxCM(t) = ⟨RxCM(t)RxCM(0)⟩/⟨(RxCM)2⟩ and ϕRzCM RzCM(t) = [⟨RzCM(t)RzCM(0)⟩ − ⟨(RzCM)2⟩]/[⟨(RzCM)2⟩ − ⟨(RzCM)2⟩], Figure 6b. Here R⃗ CM is measured for each chain

τAA ∝ N Δ , Δ ≈ 3.7 4386

(19)

dx.doi.org/10.1021/ma202745b | Macromolecules 2012, 45, 4381−4393

Macromolecules

Article

Figure 7. (a) Log−log plot of chain relaxation times in polymer brushes versus chain length N, as obtained from estimates of the asymptotic exponential decay in the MD simulations, for A = Rez2, Rex2, Rgz2, Rgxy2. Note that Rex2 is the largest relaxation time; all other relaxation times coincide within their statistical errors. The upper straight lines illustrate a power law τAA∝NΔ with an effective exponent Δ = 3.7. All data are from the MD simulation only. (b) Same as part a, but including now also the initial relaxation times for the xy and z components of the gyration radius, as well as relaxation times extracted from monomer mean square displacements (see below). Again straight lines and their slopes are indicated. (c) Initial relaxation times of the gyration radius relaxation, plotted in log−log format against N, as resulting from the Monte Carlo simulation. The right inset shows a log−log plot of ⟨Rgz2⟩ and of the fluctuation ⟨Rgz2⟩ − ⟨Rgz⟩2 = δ(Rgz2) versus N. The upper left inset shows the corresponding 2 “effective friction coefficient” ζef f = τ(i) z /Nδ(Rgz ). (d) Log−log plot of relaxation times normalized by the product of chain length N and associated static fluctuation, as obtained from MD simulations.

and for the z component and the xy components of the radii the same exponent applies, in contrast to eqs 1, 5, which imply different exponents. Remarkably, it is even the late stage relaxation times exp − τx, MSD − τx associated with xy components of the end-to-end vector which is consistently the largest relaxation time. The fact that the exponent Δ exceeds the theoretical value of eq 1, namely Δ = 3, could be attributed to a possible onset of entanglement effects (although one would expect to see strong curvature on the log−log plot in this case, which is hardly evident; but clearly relaxation times as large as 106 MD time units are extremely hard to “measure” reliably in a simulation, and possibly a crossover is not seen because our results for large N still are systematically too low). Nevertheless, the disagreement between our results for τAA and the theory14 is striking. However, Figure 7b proposes a resolution to this puzzle, including our estimates for the initial relaxation time τ(i) AA, as well as relaxation times extracted from the time evolution of monomer mean square displacements (see below) and provides thus lower bounds to τAA (these lower bounds are also consistent with eq 19). In contrast, for the initial relaxation times, we find 2 2 (i) τ(i) AA ∝ N , if A = Rgxy , in perfect agreement with eq 5, and τAA ∝ 3 2 N , if A = Rgz , in perfect agreement with eq 1. The simulation data from MC tell similar story (Figure 7c). Indeed, when we normalize the initial relaxation times τ(i) AA by the product of chain length and the associated static fluctuation (Figure 7d and the

upper left inset to Figure 7c), we should obtain an effective monomeric friction coefficient ζ, if eqs 1 and 5 are valid: (i) ξ eff ≡ τAA /[N (⟨A2 ⟩ − ⟨A⟩2 )]

(20)

One can see that for the initial relaxation times τ(i) AA the data settle down toward a plateau, so eq 20 makes sense, while with respect to the asymptotic relaxation time τAA the strong increase of ξeff with N indicates that such an interpretation is not valid. Note also that ξef f for z components and xy components in Figure 7c,d, presumably converge to almost the same value for large N. Indeed, in the monomeric friction for not too strong stretching we do not expect strong anisotropy due to the fact, that stretching does not play a role on length scales smaller than the blob size. Inside the blobs the system behaves unperturbed and can be considered isotropic. Therefore, the friction should be isotropic, too. So the theoretical concepts of Klushin and Skvortsov,14 outlined in the Introduction do work for the initial relaxation, but not for the asymptotic decay of ϕAA(t). Indeed, it is plausible that the random Brownian motions of the chains in the fluctuating potentials, exerted by their neighboring chains, will lead to a coupling of all slow degrees of freedom. So the asymptotic decay of ϕAA(t) is always controlled by the slowest relaxing variable, irrespective which property A is considered. One cannot strictly decouple the fluctuations of z component from those of x and y components, as implicitly assumed in eqs 1, 5. Possibly this behavior is due to a coupling of the chain relaxation with long wavelength 4387

dx.doi.org/10.1021/ma202745b | Macromolecules 2012, 45, 4381−4393

Macromolecules

Article

Figure 8. (a) Mean square displacement (MSD) ⟨[xi(t) − xi(0)]2⟩ of individual monomers, from i = 1 (adjacent to the grafting site) to i = 64 (chain end) for N = 64, plotted on a log−log plot as a function of time. Note that i = 36 (not shown) behaves similar to the average over all monomers of the chain. (b) Same as part a, but for ⟨[zi(t) − zi(0)]2⟩ vs time t. (c) Relaxation times τi plotted vs monomer index i, defining τi from the condition that the MSD has reached 2/3 of its saturation value. (d) Mobility of monomers parallel and perpendicular to the grafting plane.

largest fluctuation. This observation again shows that there is no simple proportionality between the relaxation time and the associated static fluctuation, as was assumed in eqs 1, 5. Of course, monomers in the middle of the chain are already rather mobile, in contrast to monomers near the grafting site, but their motions are also strongly coupled to their environment. In contrast, the end monomer which can move much more freely, experiences less coupling to its environment because most of the time it moves in a low density region. An interesting question concerning the peculiar behavior of the relaxation times τi of the individual monomers i, (i = 1, ..., N) pertains to the extent to which this behavior can be attributed to the local friction coefficient ζefi f. In order to extract these local friction coefficients, we note that the initial “ballistic regime“ (MSD ∝ t2) crosses at about t ≈ 0.5 MD time units over to a diffusive regime, where ⟨[xi(t) − xi(0)]2⟩ = 2Dxi t and ⟨[zi(t) − zi(0)]2⟩ = 2Dzi t holds in a time region 1 ≤ t ≤ 2 t.u., before the crossover to the well-known singular behavior19 MSD ∝ t1/2 [using the framework of the standard Rouse model19] sets in. These local mobilities Dxi , Dzi , extracted from a fit of the MSD in the quoted time intervals, are shown in Figure 8d. Evidently, these local mobilities are strongly reduced for a few monomers near the grafting site, and somewhat enhanced for a few monomers near the free end (the latter is more mobile, of course, since it is bound only to a single neighbor, unlike other monomers). However, for N = 64 (Figure 8(d)), one can readily verify that Dxi , Dzi are practically independent of

collective density fluctuations in the polymer brush (which are more pronounced in xy directions rather than in z directions.) Such collective fluctuations were studied by Marko et al.56,57 The fact that the concepts described in the Introduction are too much simplified is also evident when one obtains a more local view on the dynamics of monomeric motions in polymer brushes provided by the study of mean square displacements of individual monomers as a function of time (Figure 8a,b). Note that for simplicity the regime of very early times has been omitted. The mean square displacement increases with time and then saturates at a plateau value. The plateau value becomes larger, the closer the monomer index i comes to the free end; of course, since the chain is rigidly grafted, a crossover to free diffusion (as in melts) can never occur. Thus, qualitatively, this behavior is similar to mean square displacements taken in the center of mass coordinate system of a freely diffusing chain in a solution or melt.53 In the latter case, it was found convenient to define a relaxation time from the criterion that the mean square displacement (MSD) has reached 2/3 of its final saturation value.53 Of course, the number 2/3 is an arbitrary convention, so using 3/4 or 5/6, for instance, would only change the scale for such a relaxation time somewhat, but would not change general features.53 Although the plateau value increases monotonically with i, the relaxation time defined in the way just described does not (Figure 8c): The largest relaxation time occurs somewhere near the middle of the chain, and not at the free chain end although the latter clearly has the 4388

dx.doi.org/10.1021/ma202745b | Macromolecules 2012, 45, 4381−4393

Macromolecules

Article

Figure 9. (a) Log−log plot of the asymptotic relaxation times τAA versus N, comparing results for ring polymers and linear chains. Here xy and z components of Re2 are used (remember that in the case of ring polymers one considers the displacement of the monomer in the middle of the ring, i = NR/2 + 1). Part b compares autocorrelation functions for ring polymers (R) with corresponding results for linear chains (L), according to the MC simulations for early times (ϕAA(t) > 0.1).

i for 6 ≤ i ≤ 64 (note that Dxi is slightly larger that Dzi ). These inverse short time diffusivities can be identified as local friction coefficients, ζefi,xf = 1/Dxi , ζefi,zf = 1/Dzi , and thus we may conclude that the non-monotonic behavior of the local relaxation times τxi is not caused by a corresponding behavior of the local friction coefficients. Such a non-monotonic behavior is not predicted by the extension of the Rouse model to describe the dynamics of monomers in polymer brushes, using the effective potential of self-consistent field theory.17 It is also interesting that the relaxation time for the transverse displacements is larger than relaxation time for the displacements in z direction, although the magnitude of the MSDs scales in the opposite way (∝N2 in the z direction, and ∝ N in the x direction). The findings, shown in Figure 8c, clearly call for a more detailed theory of chain dynamics in polymer brushes. Finally, Figure 9a compares the asymptotic relaxation times of linear chains to the relaxation times of rings. Remarkably, all relaxation times of rings are smaller than those of the linear chains although in the initial behavior, cf. Figure 9b rings relax slower than linear chains. Additionally, the x and z components of the relaxation time of the rings do not seem to scale with the same effective exponent. While this exponent is found to be Δx ≈ 3.2 for the x component, it seems to be slightly bigger (Δz ≈ 3.5) for the z component. However, the same concerns raised above about the difficulty of measuring relaxation times of order 106 and the onset of the entanglement might be relevant here, too.

thought to be concentrated in its center of mass, that is bound to its average equilibrium position in the brush by an effective elastic spring. The strength of this elastic spring is adjusted such that the proper gyration radius components ⟨Rgz2⟩ in the z direction normal to the substrate and ⟨R2gxy⟩ in the xy directions parallel to the substrate are reproduced. In the spirit of the Rouse model, the effective friction coefficient for these center of mass motions is ζ ∝ N. Hence one predicts relaxation times that scale with different powers of N, depending on whether one considers motions perpendicular to the surface (τz ∝ N3 eq 1) or parallel to it (τxy ∝ N2, eq 5). The anomalous large exponent of τz must not be confused with the reptation law τ ∝ N3 for entangled long chains in bulk melts, rather it is due to anomalous large (“critical”14) fluctuations, ⟨Rgz2⟩ ∝ N2. These fluctuations result from the fact that chains in brushes are strongly stretched. Hence the average distance of the center of mass from the surface shows the same scaling as the brush height h, ⟨RzCM⟩ ∝ h < N, and the root-mean-square fluctuation of this distance is of the same order as the distance itself. Similarly, the free end of the chain is not always near the outer end of the brush, but rather fluctuates over the full height, 0 < Rez < h: when the free end is close to z = 0 almost all monomers must be at distances z < h/2. When the free end is close to z ≈ h, however, many monomers are at z > h/2, and hence large fluctuations in the center of mass distance from the substrate must occur. Although this prediction has been made 20 years ago, it has not yet been tested carefully to our knowledge. We have performed large scale molecular dynamics and Monte Carlo simulations, and we have found that the Klushin−Skvortsov predictions can only describe the initial relaxation of time autocorrelation functions that characterize the dynamic correlations of fluctuations of the mean square gyration radius components. We find that these autocorrelation functions show a decay with time that differs dramatically from a single exponential decay. Thus, extracting relaxation times by the frequently used method, where one asks over which time a normalized autocorrelation function has decayed from 1 to 1/e, fails completely to characterize the relaxation behavior (Figure 6) reliably. Depending on the time interval that is considered, one finds much larger effective relaxation times (Figure 6a,c) than initially. Hence the estimation of the exponential relaxation exp(−t)/τ(f)) that finally takes over as t→∞ requires a huge numerical effort. We have found that the scaling of these final relaxation times τ(f) is compatible with a power law τ(f) ∝ NΔ

V. CONCLUSIONS In the present paper we have addressed the dynamics of single polymers in polymer brushes, formed from either linear chains or ring polymers densely grafted on flat repulsive substrate surfaces under good solvent conditions. We have focused on understanding both the relaxation behavior of single polymers and their individual monomers in thermal equilibrium. The response of the polymer to nonequilibrium perturbations (such as shear deformation of the brush, or expulsion of the chain from the brush caused by cutting the bond by which the polymer is “anchored” at the substrate) is not considered at all, since such phenomena have frequently been studied in the literature already. Our study has been motivated by the simple theoretical argument of Klushin and Skvortsov,14 who argue that essentially the relaxation times associated with the motions of a polymer in a brush can be accounted for by a simple “dumbbell model”, where the mass of the whole chain is 4389

dx.doi.org/10.1021/ma202745b | Macromolecules 2012, 45, 4381−4393

Macromolecules

Article

but the ″effective″ exponent Δ ≈ 3.7 is independent of whether τ(z f) or τ(xyf) is considered. Indeed, in the late stages, due to dynamic couplings among the relaxing variables, it is plausible that the slowest relaxation dominates all autocorrelation functions in the same way. Since “effective” exponents describing the scaling with N, always depend somewhat on the range of N that is fitted (in other words, there exists slight curvature on the log−log plots of τ(i) vs N and τ( f) vs N), it is not clear whether or not the asymptotic regime could be reached. In fact, one can argue that, due to entanglement effects, a mechanism analogueuous to the “arm retraction mechanism” of star polymers should dominate ultimately, leading to an exponential variation of relaxation times with N. It is an open problem to ascertain at which N in polymer brushes such entanglement effects start to dominate. Interestingly, we have found that (noncatenated) ring polymers relax somewhat faster than linear ones, at late times (Figure 9a), although at earlier times the opposite behavior is found (Figure 9b). The theoretical understanding of ring polymer dynamics is even more rudimentary, which is not surprising in view of the fact that already the static properties (Figure 1,Figure 5) are less well understood in this case. Finally, an intriguing aspect is the dynamics of individual monomers (Figure 8). While the static fluctuations are larger the further away a monomer label is from the grafted one (Figure 2), the relaxation time of individual monomers has a nonmonotonic variation with the monomer position along the chain (Figure 8a). This nonmonotonic variation of the relaxation time is related to mean square displacements in x or y directions (parallel to the grafting surface), and it only shows up at late times. Extracting a mobility from the MSDs at short times, one finds that the mobility is essentially constant along the chain, apart from a few monomers close to the grafting point and the free chain end. Again we conclude that it is not straightforward to link relaxation times with static fluctuations for polymer brushes. Clearly it would be interesting to study the problem of chain dynamics in brushes varying also other parameters (grafting density, solvent quality, interaction with the substrate, curvature of the substrate, etc.). Also the inclusion of explicit solvent molecules may be worthwhile. Finally, it would be very interesting to have corresponding experimental studies (in principle this is possible by NMR methods when chains containing C13 rather than C12 are grafted, or by linking fluorescent markers to the chains, etc.). We hope that our study is a motivation for such efforts in the future.

Figure 10. Schematic description of a polymer chain with one end (i = 1) grafted at the point x = 0, y = 0 of the plane z = 0. N beads connected by harmonic bonds of length S move under the action of random forces and frictional forces (discrete Rouse model or Verdier− Stockmayer model, respectively). The open square shows the center of mass of the chain. Excluded volume and chain noncrossability are neglected (“phantom chain” description).

moment as penetrable, so that it does not have any physical effect. Then at times much shorter than the Rouse time τ1 [eq 4.2519], ζ being a friction coefficient τ1 =

gCM (t ) = ⟨[R⃗ CM(t ) − R⃗ CM(0)]2 ⟩ = 6DN t = 6(kBT /Nζ )t ,

t < τ1

(23)

At a time τMon  ζS /(3π kBT ) where a monomer has moved over a distance δ⃗ comparable to the bond length S , the center of mass has only moved over a squared distance given by 2

2

gCM (τMon) = (2/π 2)S2/N

(24)

When we also consider correlation functions ϕR R (t ) = ⟨R⃗e(t ) ·R⃗e(0)⟩/⟨Re 2⟩ e e

ϕR

APPENDIX: RELAXATION FUNCTIONS FOR THE END-TO-END VECTOR AND THE CENTER OF MASS MOTION INTERPRETED IN TERMS OF THE ROUSE MODEL In our work we consider a polymer chain with N beads (i = 1, ..., N) connected by bonds of length S , the first bead being grafted to the plane z = 0, Figure 10. We take the position of the grafted monomer r1⃗ = 0 as the coordinate origin, and are interested in the dynamics of the end-to-end distance R⃗ e and center of mass position R⃗ CM,

CMR CM

(t ) = [⟨R⃗ CM(t ) ·R⃗ CM ⟩/⟨R CM 2⟩

(25) (26)

it is clear that the initial decay of ϕReRe(t) will be faster (by a factor of order N) than the initial decay of ϕRCMRCM(t), since the mobility of a monomer is the order of N larger than the mobility of the whole chain. The correlator ϕReRe(t) is traditionally considered in the textbooks and in the framework of the discrete Rouse model, where the normal coordinates of the chain are X ⃗ p (t ) =

N i=1

(22)

the center of mass moves diffusively, the constraint of grafting an end at the origin then is not yet felt by the center of mass motion



R⃗e = rN⃗ − r1⃗ , R⃗ CM = (1/N ) ∑ ri ⃗

ζN 2 S2 3π 2kBT

1 N

N i=1

p = 0, 1, ..., N − 1 (21)

⎡ pπ (i − 1/2) ⎤ ⎥, ⎣ ⎦ N

∑ ri (⃗ t )cos⎢

(27)

and hence their correlators show simple exponential decay, ⟨X⃗ p(t)X⃗ p(0)⟩ = 3δ pq(kBT/κ p) exp(−t/τ p), with 58 κ p = 24(NkBT/l2) sin2(pπ/2(N − 1)), τp = (2Nζ/κpp2)) which yields

As a reference model, we consider the Rouse model of harmonic phantom chains, and take the surface at z = 0 for the 4390

dx.doi.org/10.1021/ma202745b | Macromolecules 2012, 45, 4381−4393

Macromolecules

Article

(for pπ ≪ 2(N − 1) one has κp ≈ 6kBT(p2π2)/l2N), simply τp = τ1/p2, τ1 = ζl2N2/(3kBTπ2). Expressing then the coordinates ri⃗ (t) back in normal coordinates, and subtracting the center-of-mass motions, the correlation function of P⃗  R⃗ e(t) − R⃗ CM can be written as ⟨P⃗ (t)P⃗ (0)⟩ = 16∑ p =1, 3,... ⟨X⃗ p (t)X⃗ p (0)⟩ = 16∑p=1,3,...(3kBT/κp) exp(−(t)/(τp)) so that ϕ R R (t ) = e e

ϕR

2const S2N[1 − ϕR

Note that in the continuum approximation to the Rouse model where differences ri⃗ +1(t) + ri⃗ −1(t) − 2ri⃗ (t) are approximated as ∂2ri⃗ (t)/∂t2, eq 28 holds as a strict equality, and the sum over p extends to ∞. For times t→∞ eq 28 is dominated by a single exponential: ϕR R (t ) ≡ (8/π 2) exp(−t /τ1)

e e

t (8/π 2) τ1

⎛ 4N ⎞ 1 ≈ 1 − t⎜ ⎟ ⎝ π Ωτ1 ⎠ p = 1,3,...



(30)

Note that at this point it is crucial that N is finite and the number of normal modes [eq 27] also is finite. In order to understand the limiting behavior as N → ∞, it is convenient to consider the time derivative of ϕReRe(t), d ϕR R (t ) = (8/τ1π 2) dt e e

t ≪ τ1

ϕ R zR z(t ) = [⟨Rez(t )Rez(0)⟩ − ⟨Rez⟩2 ]/[⟨(R ez)2 ⟩ − ⟨R ez⟩] e

e

(37)



2

exp(−tp /τ1)

ϕR z

Z CMR CM

(31)

p = 1,3

∫0



(32)

so in this limit N → ∞ ϕReRe(t) develops a square root behavior and the initial relaxation time cannot be defined ϕR R (t ) ≈ 1 − (4/π 3/2)(t /τ1)1/2 , e e

N→∞

zz z z gCM (t ) =⟨[R CM (t ) − R CM (0)]2 ⟩ = (2kBT /Nζeffz )t , eff t < τ1

t /τ1 ≪ 1, (33)

xyxy x x gCM (t ) =⟨[R CM (t ) − ⟨R CM (0)⟩]2 y y + ⟨[R CM (t ) − R CM (0)]2 ⟩

For finite N the behavior is analytic at short times, since eq 28 is a finite sum of exponentials, and hence analytic at t = 0. But the crossover from eq 30 to eq 33 already occurs at times t of the order of tcross = (πτ1)/N2 = ζS 2/(3πkBT) = πτMon. So the (i) initial relaxation time τ(i) ReRe defined from ϕReRe(t) = 1 − t/τReRe as t→0, τR(ie)R e = τ1π 2/(4N ) = 3N S2/(12kBT )

(38)

Although now eqs 30, 33, and 35 are no longer quantitatively accurate, in qualitative respects the previous conclusions remain unchanged. (i) The conclusion that τRzCM RzCM must scale like N3 rather than N2 in a polymer brush can then be obtained by generalizing the above treatment as follows: 2 (I) eq 23 is replaced by (τeff 1 is a time presumably ∝ N , as for a free chain)

dp exp( −tp2 /τ1)

= (2/π 3/2)(tτ1)−1/2

z z z 2 (t ) = [⟨R CM (t )R CM (0)⟩ − ⟨R CM ⟩] z 2 z 2 ⟩] /[⟨(R CM ) ⟩ − ⟨R CM

for large N and small t the sum over the Rouse modes can be approximated by an integral d ϕ (t ) ≈ (4/τ1π 2) dt R eR e

(t )] = 6(kBT /Nζ )t ,

which is equivalent to eq 35. Of course, when one makes the assumption (closer to the situation considered in the simulation) that the surface is impenetrable rather then penetrable as assumed so far, an anisotropy between z components and xy components is introduced into the system. Since now the vectors R⃗ e, R⃗ CM can have only positive z components, we can no longer conclude that ⟨R⃗ e⟩ = 0, ⟨R⃗ CM⟩ = 0, as implicitly assumed in defining the functions ϕReRe(t), ϕRCMRCM(t) in eqs 25, 26. Such a symmetry with respect to the coordinate origin (the grafting site) of the chain still holds with respect to the x and y components, for which correlators are defined as in eqs 25 and 26, while for the z components the proper definitions are

while for small t we have (there are about N/2 “odd Rouse modes” p = 1, 3, ...) ϕ R R (t ) ≈ 1 −

CMR CM

(36)

(29)

e e

(35)

i.e., in the Rouse model the initial relaxation time of the center of mass correlator is proportional to the Rouse time itself. We note that eq 35 can be trivially justified by taking the square in eq 23 and noting that for our choice of coordinate system ⟨[R⃗ CM(t)]2⟩ = ⟨[R⃗ CM(0)]2⟩ = const S2N and hence

(28)

p = 1,3, ···

) (t ) = 1 − t /τR(iCM R CM ,

) τR(iCM ∝ ζN 2 S2/kBT ∝ τ1 ,

(8/p2 π 2) exp(−tp2 /τ1)



CMR CM

=(4kBT / Nζeffxy)t ,

(39)

t < τ1eff

τeff 1

(40)

Thus, we assume that a time ≫ = 2 xy τxy Mon = ζeff/(3π kBT) exists, up to which the center of mass motion still can be described as a Rouse-like diffusion (possibly friction coefficients ζzeff in z direction differ from the friction coefficient in xy directions and from the friction coefficient for a free chain (ζ)). Now eqs 39 and 40 can be rearranged in analogy with eq 36 as follows

(34)

rules the behavior of ϕReRe(t) not up to t ≈ τ(i) ReRe but only for times of the order of the “microscopic” time τMon. Thus, for ϕ(i) ReRe the concept of the initial relaxation time is not very useful in practice. However, with respect to the correlator describing the center of mass motion, eq 25, the situation if clearly different, since we have found above that the initial decay of ϕRCMRCM(t) with time t is about a factor of N smaller than the initial decay of ϕReRe(t). As a consequence, we suggest that

τzMon

ζzeff/(3π2kBT),

z 2 z z 2[⟨R CM ) ⟩ − ⟨R CM (t )R CM (0)⟩] = (2kBT /Nζeffz )t z 2 z 2 2[⟨R CM ) ⟩ − ⟨R CM ⟩ ][1 − ϕ R z

z CMR CM

(41)

(t )] = (2kBT /Nζeffz )t (42)

[⟨(RzCM)2⟩

Using then that − ∝ N in a strongly stretched brush, we have made it plausible that 4391

(RzCM)2]

2

dx.doi.org/10.1021/ma202745b | Macromolecules 2012, 45, 4381−4393

Macromolecules ϕR z

z CMR CM

Article

(t ) = 1 − t /τ R(iZ)

z CMR CM

t < τ1eff

,

(13) Spirin, L.; Galuschko, A.; Kreer, T.; Binder, K.; Baschnagel, J. Phys. Rev. Lett. 2011, 106, 168301. (14) Klushin, L. I.; Skvortsov, A. M. Macromlecules 1991, 24, 1549. (15) Binder, K.; Lai, P.-Y.; Wittmer, J. Faraday Discuss. 1994, 98, 97. (16) Neelov, I. M.; Binder, K. Macromol. Theory Simul. 1995, 4, 1063. (17) Johner, A.; Joanny, J.-F. J. Chem. Phys. 1993, 98, 1647. (18) He, G. L.; Merlitz, H.; Sommer, J.-U.; Wu, C.-X. Macromolecules 2007, 40, 6521. (19) Doi, M.; Edwards, S. F. The Theory of Polymer Dynamics; Clarendon Press; Oxford, U.K. , 1986. (20) Kuhn, W.; Kuhn, H. Helv. Chim. Acta 1945, 28, 1533. ibid. 1946, 29, 71. (21) Cosgrove, T.; Heath, T.; van Lent, B.; Leermakers, F.; Scheutjens, J. M. H. S. Macromolecules 1987, 20, 1692. (22) Milner, S. T.; Witten, T. A.; Cates, M. E. Macromolecules 1988, 21, 2610. (23) Skvortsov, A. M.; Pavlushkov, I. V.; Zhulina, Y. B.; Borisov, O. V.; Pryamitsyn, V. A. Polym. Sci. U.S.S.R. 1988, 30, 1706. (24) Murat, M.; Grest, G. S. Phys. Rev. Lett. 1989, 63, 1074. (25) Murat, M.; Grest, G. S. Macromolecules 1989, 22, 4054. (26) Lai, P.-Y.; Binder, K. J. Chem. Phys. 1991, 95, 9288. (27) Wittmer, J.; Johner, A.; Joanny, J.-F.; Binder, K. J. Chem. Phys. 1994, 101, 4379. (28) Kreer, T.; Metzger, S.; Müller, M.; Binder, K.; Baschnagel, J. J. Chem. Phys. 2004, 120, 4012. (29) Dimitrov, D. I.; Milchev, A.; Binder, K. J. Chem. Phys. 2007, 127, 084905. (30) DeGennes, P. G. Scaling Concepts in Polymer Physics; Cornell Univ. Press: Ithaca, NY, 1979. (31) Reith, D.; Milchev, A.; Virnau, P.; Binder, K. EPL 2011, 95, 28003. (32) Weeks, J. D.; Chandler, D.; Andersen, H. C. J. Chem. Phys. 1971, 54, 5237. (33) Kremer, K.; Grest, G. S. J. Chem. Phys. 1990, 92, 5057. (34) Milchev, A.; Paul, W.; Binder, K. J. Chem. Phys. 1993, 99, 4786. (35) Milchev, A.; Binder, K. Macromol. Theory Simul. 1994, 3, 915. (36) Milchev, A.; Binder, K. J. Comput.-Aided Mater. Des. 2002, 9, 33. (37) Binder, K.; Paul, W. J. Polym. Sci., Part B: Polym. Phys. 1997, 35, 1. (38) Anderson, J.; Lorenz, C.; Travesset, A. J. Comput. Phys. 2008, 227, 5342. (39) Soddemann, T.; Dünweg, B.; Kremer, K. Phys. Rev. E 2003, 68, 046702. (40) Dünweg, B.; Stevens, M.; Kremer, K. In Monte Carlo and Molecular, Dynamics Simulations in Polymer Science; Binder, K., Ed.; Oxford Univ. Press: New York, 1995; p 125. (41) Yin, F.; Bedrov, D.; Smith, G. D.; Kilbey, S. M. J. Chem. Phys. 2007, 127, 084910. (42) Cates, M. E.; Deutsch, J. J. Phys. (Fr.) 1986, 47, 2121. (43) Müller, M.; Wittmer, J. P.; Cates, M. E. Phys. Rev. E 1996, 53, 5063. (44) Brown, S.; Szamel, G. J. Chem. Phys. 1998, 108, 4705. ibid. 109, 6184. (45) Müller, M.; Wittmer, J. P.; Barrat, J. L. Europhys. Lett. 2000, 52, 406. (46) Müller, M.; Wittmer, J. P.; Cates, M. E. Phys. Rev. E 2000, 61, 4078. (47) Hur, K.; Winkler, R. G.; Yoon, D. Y. Macromolecules 2006, 39, 3975. (48) Vettorel, T.; Grosberg, A. Y.; Kremer, K. Phys. Biol. 2009, 6, 025013. (49) Hur, K.; Jeong, C.; Winkler, R. G.; Lacevic, N.; Gel, R. H.; Yoon, D. Y. Macromolecules 2011, 44, 2311. (50) Halverson, J. D.; Lee, W. B.; Grest, G. S.; Grosberg, A. Y.; Kremer, K. J. Chem. Phys. 2011, 134, 204904−204905. (51) Reith, D.; Mirny, L. A.; Virnau, P. Prog. Theor. Phys. Suppl. 2011, 191, 135. (52) Paul, W.; Binder, K.; Heermann, D. W.; Kremer, K. J. Phys. II (Fr.) 1991, 1, 37.

(43)

with z 3 2 eff ) τR(iCM z z ∝ ζeff N S / kBT (≫τ1 ) R CM

(44)

It is straightforward to show analogously from eq 40 that ϕ R xy R xy (t ) = 1 − t /τR(i)xy R xy , CM CM

CM CM

t < τ1eff

(45)

with τR(i)xy R xy ∝ ζeffxyN 2 S2/kBT CM CM

(46)

However, these simple arguments do not tell what controls the time scale τeff 1 over which free Rouse-like motion of the center of mass can occur: (i) Coupling between longitudinal and transverse components of the center of mass displacements? (ii) Cooperative motions of several chains? (iii) Entanglement effects, etc.? Of course, it would be desirable to carry out a more complete calculation studying the Rouse model for a polymer chain in a brush, where the presence of the other chains is described self-consistently in terms of a potential U(z) = kBTπ2(h2 − z2)/(8S2N 2), h being the height of the brush, following the approach by Johner and Joanny.17 These authors have shown that the asymptotic relaxation time τRzCM RzCM scales as N3 as well, in the framework of the Rouse description. However, a detailed (numerical) analysis of such a model must be left to future work.



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work is supported by Deutsche Forschungsgemeinschaft (DFG), grant numbers SFB625 A3 and A17, and Bi 314/23. Computing time on the GPU cluster of the Center for Computational Sciences Mainz at the ZDV Mainz is gratefully acknowledged.



REFERENCES

(1) Halperin, A.; Tirrell, M.; Lodge, T. P. Adv. Polym. Sci. 1992, 100, 31. (2) Grest, G. S.; Murat, M. In Monte Carlo and Molecular Dynamics Simulations in Polymer Science; Binder, K., Ed.; Oxford Univ. Press: New York, 1995; p 476. (3) Klein, J. Annu. Rev. Mater. Sci. 1996, 26, 581. (4) Szleifer, I.; Carignano, M. A. Adv. Chem. Phys. 1996, 94, 165. (5) Grest, G. S. Adv. Polym. Sci. 1999, 138, 149. (6) Leger, L.; Raphael, E.; Hervet, H. Adv. Polym. Sci. 1999, 138, 185. (7) Advincula, R. C., Brittain, W. J., Caster, K. C., Ruehe, J., Eds. Polymer Brushes,; Wiley-VCH: Weinheim, Germany, 2004. (8) Descas, R.; Sommer, J.-U.; Blumen, A. Macromol. Theory Simul. 2008, 17, 429. (9) Binder, K.; Kreer, T.; Milchev, A. Soft Matter 2011, 7, 7159. (10) Egorov, S. A.; Milchev, A.; Klushin, L.; Binder, K. Soft Matter 2011, 7, 5669. (11) Galuschko, A.; Spirin, L.; Kreer, T.; Johner, A.; Pastorino, C.; Wittmer, J.; Baschnagel, J. Langmuir 2010, 26, 6418. (12) Spirin, L.; Galuschko, A.; Kreer, T.; Johner, A.; Baschnagel, J.; Binder, K. Eur. Phys. J.E 2010, 33, 307. 4392

dx.doi.org/10.1021/ma202745b | Macromolecules 2012, 45, 4381−4393

Macromolecules

Article

(53) Paul, W.; Binder, K.; Herrmann, D. W.; Kremer, K. J. Chem. Phys. 1991, 95, 7726. (54) Alexander, S. J. Phys. (Fr.) 1977, 38, 977. (55) De Gennes, P. G. Macromolecules 1980, 13, 1069. (56) Marko, J. F.; Witten, T. A. Macromolecules 1992, 25, 296. (57) Marko, J. F.; Chakrabarti, A. Phys. Rev. E 1993, 48, 2739. (58) Verdier, P. H.; Stockmayer, W. H. J. Chem. Phys. 1962, 36, 227.

4393

dx.doi.org/10.1021/ma202745b | Macromolecules 2012, 45, 4381−4393