Faraday Symp. Chem. SOC.,1983 18 37-47 Dynamics of Entangled Flexible Polymers Monte Carlo Simulations and their Interpretation BY ARTUR BAUMGARTNER AND KURT BINDER* KURTKREMER Institut fur Festkorperforschung der Kernforschungsanlage Julich D-5170 Julich Postfach 1913 West Germany Received 1st August 1983 Monte Carlo simulations have been reviewed for several models of dense polymer systems (i) a system of very short (N = 16 links) freely joined chains with Lennard-Jones interactions in the continuum (ii) self-avoiding walks on a diamond lattice with up to N = 200 links and (iii) pearl-necklace chains with up to N = 98 hard spheres. The cases of both a mobile and a frozen-in environment have been considered. The dynamics of displacements and the structure factors have been studied and interpreted in terms of the Rouse model and the reptation model paying particular attention to the crossover between these models.Only for the frozen- environment case is reptation fully verified while for models (i) and (ii) simple Rouse behaviour is found and for model (iii) relaxation with a diffusion constant D cc N-2-o*o.2and disengage- ment time zdcc 1v3.4*0.4is found but the monomer displacements are in disagreement with reptation laws. A brief comparison with pertinent experiments has been made and the crossover between the Rouse model and the Zimm model is considered briefly. 1. INTRODUCTION CURRENT MODELS FOR CHAIN DYNAMICS Static properties of long flexible polymers in dilute and semi-dilute solutions are understood by renormalization group methods,l and in concentrated solutions or melts they follow trivial random-walk statistics (radius R x IN; for chains of N links of length 1).2 Dynamic properties in contrast are less well understood and current theoretical discussions rest on very simplified model~.~-l~ In the Rouse the forces of the environment on a chain are represented by a heat bath inducing local conformational changes.The characteristic time zq over which fluctuations characterized by a wavevector q decay is then3? 7* where W is the rate at which links change their orientation and v x 0.59 for chains which are swollen due to excluded-volume interactions.' Eqn (1 a) describes the ~ relaxation for times t shorter than the chain-relaxation time T~ oc/ W-1N2v+1 while /~ for t > T~ diffusive behaviour eqn (1 b) dominates the motion for qR > 1.20 Similarly the monomer displacements (r%(t))behave as3 37 DYNAMICS OF ENTANGLED POLYMERS In the Zimm the hydrodynamic backflow forces are also included.This is important for dilute solutions where W cc qol,qo being the solvent viscosity and Instead of eqn (2a)for the monomer displacements we now have (rt(t)) x P(Wt); for t < z,/~. These displacements show up in the 'incoherent scattering factor' &,(q t) cc (1/N)Z (exp {h [(ri(t)-ri(O)ll<% ~XP [-+q2(ra(t))l* For both the Rouse and the Zimm model the coherent structure factor Sco,(q 0 cc N-2c (exp {iq [ri(t)-q(0)lD (54 ij takes the simple form &oh(q t) cc exp [-const (qZ)2(Wt).] Wt 9 1 (5b) where the exponent n = 1 in the diffusive regime while for qR > 1 and t < z,,~ one has n = (Rouse model) or n = (Zimm model).The 'reptation model' of de Gennes97 l23 l5 and EdwardslO? l3considers dense systems where the chains are entangled. It is thought that the effect of these entanglements is to create for each chain a tube of diameter dT,such that the chains are restricted to snake-like motions along the tube axis. On the local scale inside the tube i.e. for (ra(t)) 5 d& the chain performs Rouse-like motions [eqn (2)].For larger times one has insteadg idT( Wt)' W-'(dT/l)* < t <ZRouse OC W-lP (64 (rt(t)) % /dT(Wt/N)i t~~~~~ < t < zd cc (Z/~T)~N~ (6 b) [Dt t > zd D cc Wd+NP2. (64 While this behaviour also shows up in Si,,(q t) [eqn (4)] &h(q t) is more complicated as there is no single time zq dominating the scattering under wavevector q.15For R-l < q < G1it is predicted that Scoh (q,t)decays non-exponentially towards a constant value15 Scoh (4,t)/scoh (q' O) -q2 d+/36 < zd (7a) while later the relaxation should be indep endent of q Scoh (4,t)/scoh (4,O) exp (-t/zd ) > zd.(7b) All these models oversimplify the actual dynamics which are those of a dense fluid in the local environment. It is not completely clear under which conditions these models are valid. In addition one needs a better microscopic understanding of parameters such as W and dT.In fact while reptation seems to be a nice framework for bulk viscoelastic properties,13 the precise microscopic meaning of the entanglements is not so clear.We would also like to study the crossover between the various models which occurs when N or the concentration c of a polymer solution is varied.11y12q21 Even the derivation of some of the results is rather qualitative and doubts can be raised about their validity; e.g. it has been suggested20 that instead of eqn (7b) one should have s,,h(q t)/s,,h (q,O) cc exp (-Dq2 t) for t 9 zd as in the case qR < 1. Since Monte Carlo simulations have been useful for understanding static chain properties and the asymptotic laws valid for N +co could be verified for short chains,22 one expects simulations to be useful for chain dynamics also. Such A. BAUMGARTNER K.KREMER AND K. BINDER 39 simulations have indeed been performed ~ecently~~-~~ and have a bearing on the above models as well as on neutron-scattering experiment^.^^^ 28 In the following we have summarized these studies and pointed out some questions which still need answers. 2. SIMULATIONS OF SHORT CHAINS IN THE CONTINUUM A simple model consists of rigid links of length 1 freely jointed together at arbitrary angles. Interactions can be introduced for example by postulating a Lennard-Jones potential U(rij)= 4&[(a/rij)l2-(c~/r~~)~] between any pairs of beads at points ri and ri [rij= ri-ri3. While E determines the temperature scale a= 0.41 was chosen23 since then the static properties follow asymptotic laws down to very small values of N.29 Dynamics are introduced by randomly choosing a bead of chain i and moving it through a randomly chosen angle 4 on a circle while keeping all other bead positions fixed.This trial move is accepted only if the transition probability W(ri-r;)exceeds a random number q-with 0 < q-< 1 otherwise it is rejected and another move is tried. A transition probability which both satisfies detailed balance with the equilibrium probability distribution Po cc exp (-Z/kT) where X = ZU(rij)is the Hamiltonian of the system and simulates entanglement restrictions is [SZ = Wj # i ri) -Z(rj # i 4)l exp(SX/kT) if dX' < 0 no intersection if &%>O (0 if the move would require link intersection. (8 c> In our simulation of melts eqn (8 c) reduces the rate of accepted moves to ca.one-fifth of what it would be with eqn (8a)or (8b)alone. In the related work of Bishop et ~1.~~ longer chains [but only eqn (8a)and (8b)lare used; otherwise their results are similar to those of ref. (23). As a first step this model was studied in the dilute limit of isolated single chains with N = 16 links. It was found that eqn (l) (2) (4) and (5) account for the data. Only for very large q (41 1)does the gaussian approximation for ri(t)-r,(O)involved in eqn (4) become invalid as expected. As a second step 10 such chains were put in a box of size L = 41 and periodic boundary conditions were applied to simulate a macroscopic system. Again the results are in agreement with the Rouse model cf. fig. 1. Displacements (rf(t)) and SinC(q, t) show distinct the Rouse behaviour [eqn (2a)l.No intermediate reptation regime [eqn (6a)l is seen. Also Sc,,(q t)is in quantitative accord with Rouse behaviour. At longer times one sees crossover to diffusion of the chains as a whole. (Note the difference between Scohin the laboratory system and in the centre-of-gravity system in fig. 1; the crossover occurs at t z zIlRindependent of q in accord with theory.20) As a third step after equilibrating this system the configurations of all but one of the chains were frozen-in and only one chain allowed to move. Thus we simulated the situation of one chain moving in the presence of randomly fixed obstacle^.^ Now the predicted behaviour (rg(t)) cc ti [eqn (6a)lwas readily seen over at least three decades of time and in the same time interval the centre of gravity rcgfollowed the law9(r&(t)) cc ti.These results prove that in spite of the shortness of the chains there must be many entanglements along each chain. Furthermore the geometry of these entanglements at the time when we stop the movements of all other chains does not change substantially. Hence the fact that in the frozen-in case we clearly see reptation and in the mobile case we do not show that dT should not be linked exclusively to 40 DYNAMICS OF ENTANGLED POLYMERS 30 2 10 1.3 1 0.7 - n k W 01 a I& .5 0, Y c ‘i 1 lo-’ 0.7 1 0.L -2 v c \ 216’ n 4 G W c0 6 5 0 0 1 1o3 lo4 * 1o5 lo6 Fig.1. Log-log plot of S,,,(q t) (upper part) and &,h(q f)/&& 0) (lower part) against time (in units of attempted moves per bead) for various q. Solid curves represent the Rouse model prediction where the rate W[eqn (l)] was adjusted (W x 0.025 in our time units). Open circles represent the structure factor in a coordinate system in which the centre of gravity of the considered chain would be at rest. c = 2.5 kT/&= 3. From ref. (23). geometrical constraints of chain configurations but should also take into account the local chain mobility. Of course one may quote many reasons why melt simulations fail to see reptation (i) chain lengths N too short (ii) chain density too low (iii) temperature too high or (iv) reptation model inadequate.In order to check for (ii) and (iii) simulations were also performed at a four times higher density and at kT/c =0.4 respectively. In both cases the chains were found to be frozen into a sort of glass-like state with displacements (rf(t)) < 12 over the timescales of interest. 3. SIMULATIONS OF LONG LATTICE CHAINS Self-avoiding walks (SAW) on the diamond lattice with up to N =200 steps were simulated2s at a concentration of c x 0.344. The analysis of static properties revealed that on a small scale (within blobs of size NBx 20) the chains are still ‘swollen’ while on a larger scale the chains behave like ideal random walks. However simple Rouse A. BAUMGARTNER K. KREMER AND K. BINDER t h' 1oo l"l(l r IIIIIIIII loo lo2 104 lo6 t Fig.2. Log-log plot of gl(t) = (1/20) Z& (rf(t)) against time the latter measured in attempted moves per bond for chains with N = 200 at the tetrahedral lattice and a concentration c = 0.344.g2(r)is the same mean-square displacement of inner monomers of the chain measured in the centre-of-gravity systems of each chain while g3(r)= (r&(t)).From ref. (25). 1.0 .op, I 0.o 10.0 Fig. 3. Coherent structure factor s,,h(q t)/&,h(q 0) for a system of chains with 200 links moving on a tetrahedral lattice at a concentration c = 0.344. Data are plotted in semilog form against the variable (1/6)q21 2/ Wr to verify the asymptotic decay proportional to exp (-const q2I22/Wr)predicted by the Rouse model. The parameter W is set equal to unity.k = 0,0.1; 0.3; 0 for poly(dimethy1 V 0.2; + 0.25; 0 0.4; x ,0.5; A,0.6. Experimental data (0) siloxane) melts2' and the simulation data of fig. 1 (+)are also included fitting W = 4.3 for both. From ref. (25). behaviour is again observed even for the longest chains (fig. 2). It is seen that eqn (2a) holds until the mean-square displacement becomes of the order of the square of the gyration radius (RL) of the chain; then diffusion of the chain as a whole takes over. Also S,,,(q t) (fig. 3) confirms the Rouse-like behaviour. The simulation shows that in the regime where one probes the internal motions of the chains the normalized structure factor depends on qand t asq22/tonly. Deviations from this behaviour which lead to quicker decay are seen at long times and are due to the diffusive motion of the chains.By contrast according to the reptation model one would expect a DYNAMICS OF ENTANGLED POLYMERS 0 .A n + 4 .y v t. 0.5 t 0.01lIIrIll111IrlIIIlIlIrlll1 0.o 0.5 1.O 1.5 2.0 2.5 (&k212fh)-I Fig. 4.Coherent structure factor &,h(q t)/&,h(q 0) plotted against (&q212dWt)-l for the case of a single mobile chain in a frozen-in environment. Intercepts are used to estimate the tube diameter dT from eqn (74. From ref. (25). k S(k,t =a) dT N 0.3 0.50f0.05 275 45f 10 A 0.4 0.33f0.05 250 41 f5 + 0.5 0.21k0.02 225 38f5 0 0.6 0.13f0.02 205 34f5 l'"'~'~'~~'"'1 0 50 100 150 Fig. 5.Planar projections of 'snapshots' of chain configurations on the tetrahedral lattice confined to straight tubes of various diameters dT:(a) 32 (b) 24 and (c) 16 (N = 400,1= d3).From ref. (31). A. BAUMGARTNER K. KREMER AND K. BINDER 43 crossover to slower decay and the data should not scale as there is an additional dependence on the parameter 4dTe2' Such behaviour occurs if the environment of the moving chain is frozen-in (fig. 4). The resulting tube diameter dT =10 lattice spacings is of the same order as the length over which the SAW interaction is screened out in rough agreement with the expectation12 that there should be only one characteristic length in the system. Obviously this is not true in melts the screening length is of the order of a few A (monomer distances). Estimates for dTextracted from viscoelastic propertiesl3? 30 are dT =30-80 A.If dT were so large for real systems the experimental data27 included in fig. 3 would satisfy qdT 2 2 and might then be affected only weakly by the tube constraints. Thus it is necessary to study in more detail the crossover from Rouse to reptation behaviour. As a first step simulations of chains in straight tubes of various diameters have begun (fig. 5).31 7 4 3.0 I111111 I I11l11ll I I,11111 1 IIIIIII 1o3 104 105 106 107 Fig. 6. Log-log plot of gl(t) against time for the single mobile chain (N =200) in a frozen-in environment built up by similar chains at c =0.344. g2(t)is the mean-square displacement measured in the centre-of-gravity system. From ref. (25). While there are still uncertainties about coherent scattering from one reptating chain,l592of 239 25 the predictionss [eqn (6a)and (6b)l for the displacements are verified when one considers a chain in a frozen-in environment (fig.6):Rouse behaviour occurs for (ri(t)) 5 d+,while (rt(t)) Kti holds until (rt(t)) reaches (&),where one again finds (rp(t)) K&until at zd diffusive behaviour sets in [eqn (6c)l. Estimates for dT extracted from the various regimes are consistent with each One may again suspect of course that the failure to see reptation in the case where all chains are mobile might be due to too short chains and/or too small c~ncentration.~~ This is not the case for the experimental data in fig. 3 the molecular weight used (M =60000) is safely in the regime where the viscosity behaves as q KMk4 (this regime starts at M =1500027) and hence the chains are long enough to be in the strongly entangled regime.DYNAMICS OF ENTANGLED POLYMERS lo2 -I 1 I ,I,![ I I I I I I,lI I I I I I III 10 7 c w-b o-1 1 I I I1 1111 I I I I Illll I 1 I I1 I11 1o3 1o4 t 1o5 106 Fig. 7. Log-log plot of g(t) against time for a system of 24 chains of 'pearl-necklace' chains with N = 72 hard spheres each. Intersections of the straight lines serve to define various crossover times 7, 7; and 7d. From ref. (26). 10 N 100 I I I I IIIll1 / II 10 n c W bo" 4I'. .. . . .... .. .. .. .. . ..,,. *..* I I 1 10 102 1o3 1o4 1 t Fig. 8. Log-log plot of the mean-square displacement in the centre-of-gravity system plotted against time for several N (a) 72 (b) 32 and (c) 18.The chain-length dependence of the asymptotic value g!") is plotted in the upper left-hand corner. From ref. (26). 4. SIMULATIONS OF LONG CHAINS IN THE CONTINUUM The most extensive simulations to date26 concern ensembles of Np chains each consisting of N hard spheres of diameter h/Z = 0.9 freely joined together put in a box of size L chosen such that the concentration c= NpN/(L/Z)3= 0.7 for N up to 98 (and Npup to 28). In this case the diffusion constant for N 2 40 is indeed consistent with the reptation law [eqn (6c)l D cf N-2 as also found e~perimentally.~~ However the mean-square displacements do not follow the reptation predictions (fig. 7 and 8) A. BAUMGARTNER K.KREMER AND K. BINDER 0.3 ' I I I I I I I* 0.02 0.04 0.06 0.08 0.10 0.12 0.14 41A-I Fig. 9. Exponent n [defined in eqn (5b)]plotted against wavevector q for poly(dimethy1 siloxane) dissolved in deuterated benezene at concentrations c = 0.18 (triangles) and c = 0.45 (circles). From ref. (28). both g(t) and the squared displacement in the centre-of-gravity system g,(t) increase proportionally to F until t x zRouse where g,(t) saturates at gr(co)cc N. Note the contrast between g,(t) in fig. 8 and g2(t)in fig. 6 g,(t) does not show dependence on tf,in contrast with the case of the frozen-in environment. Since the saturation value satisfiesg,( co)cc N the time z where crossover to this value occurs must be identified as the relaxation time rRousecc W"2.However while in fig. 2 diffusion sets in at this time [and therefore D cc 1/N since DrRousex g,(aO) a N] in fig. 7 the displacements increase more slowly after this time and diffusion sets in at a much later time rd which varies with chain length as2' zd cc IV3.4*0*4. If this displacement behaviour in polymer melts persists for N +00 it implies that reptation theorye9l5 cannot apply. At present it cannot yet be excluded that the behaviour seen in fig. 7 and 8 is particular to the regime for crossover from Rouse to reptative behaviour. If so fig. 4 and 6-8 imply that this crossover behaviour in melts is different from the crossover in the case of the frozen-in environment. Standard theoryetl5 does not distinguish between these two cases. However various additional relaxation mechanisms might be important even if a tube model applies to melts the tube may locally expand or or the tube itself may show local Rouse-like diffusive motion together with the chain it contains; the topological constraints would be left invariant1°J3 under both these types of motion.It is not obvious to us that these mechanisms can be completely absorbed in an increase in the effective tube diameter dT. 5. CROSSOVER FROM ZIMM TO ROUSE BEHAVIOUR IN CONCENTRATED SOLUTIONS The hydrodynamic interactions important for dilute solutions4 are not easily incorporated into Monte Carlo simulations. It is believedll? 12,21 that hydrodynamic interactions are screened in semidilute and concentrated solutions and the Rouse model should then be applicable over a wide range of molecular weights.Then the simulations discussed in sections 2-4 might be compared to such systems. Experiments on poly(dimethy1 siloxane) dissolved in deuterated benzene at various concentrations have been carried out with the neutron spin-echo method.28 Analysing the scattering with eqn (5b)shows a crossover from (unscreened) Zimm relaxation (n = f)seen at large q to Rouse behaviour (n = &)at smaller q for low concentrations (fig. 9) as expected. But surprisingly at higher concentrations a second crossover is 46 DYNAMICS OF ENTANGLED POLYMERS detected from Rouse behaviour back to Zimm behaviour upon decreasing q. This phenomenon has been interpreted by incomplete screening of the hydrodynamic interaction.28 The average hydrodynamic forces are then modelled asz8 where ylo is the pure solvent viscosity ~(c)some effective solution viscosity and th(c) the hydrodynamic screening length.While the first term on the right-hand side of eqn (10) is standard,12 the last term is new. It leads to behaviour described by the Zimm model eqn (3a),in a regime R-l <q <[qo/q(c)]t;l(c) but with a reduced rate [W cc q-l(c)] screened Zimm relaxation.28 It is unclear how this hydrodynamic response over intermediate-length scales interferes with entanglement restrictions. 6. CONCLUSIONS Both neutron-scattering experiments and simulations reveal that the simple Rouse model has a wide range of applicability in understanding the time dependences of displacements and scattering intensity even in cases where entanglements affect the centre-of-gravity motion and/or viscous response of the polymer solution or melt.The crossover from Rouse behaviour to entangled behaviour in melts seems to be different from cases with a frozen-in environment; only in the latter case do simulations verify the reptation concepts. The need for more detailed theories for this crossover regime emerges since this regime also seems to be relevant for experiments. P. G. de Gennes Scaling Concepts in Polymer Physics (Cornell University Press Ithaca New York 1979). P. J. Flory Principles of Polymer Chemistry (Cornell University Press Ithaca New York 1967). P. E. Rouse J. Chem. Phys. 1953 21 1272. B. H. Zimm J. Chem. Phys. 1956 24 269.F. Bueche The Physical Properties of Polymers (Interscience New York 1962) and references therein. W. W. Graessley J. Chem. Phys. 1965 43 2696; 1967 47 1942. ’P. G. de Gennes Physics 1967 3 37. E. Dubois-Violette and P. G. de Gennes Physics 1967 3 181. P. G. de Gennes J. Chem. Phys. 1971 55 572. lo S. F. Edwards and J. M. V. Grant J. Phys. A 1973 6 1169 1 186. l1 K. F. Freed and S. F. Edwards J. Chem. Phys. 1974 61 3626. l2 P. G. de Gennes Macromolecules 1976 9 587 594. l3 M. Doi and S. F. Edwards J. Chem. Soc. Faraday Trans. 2 1978 74 1789; 1802; 1818; 1979 75 38. l4 J. D. Ferry Viscoelustic Properties of Polymers (Wiley New York 1980). l5 P. G. de Gennes J. Chem. Phys. 1980 72 4756; J. Phys. (Paris) 1981 42 735. l6 C. F. Curtiss and R. B. Bird J.Chem. Phys. 1981 74 2016 2026. l7 L. Leger and P. G. de Gennes Annu. Ret.. Phys. Chem. 1982 33 49. W. W. Graessley Adu. Polym. Sci. 1982 47 68. l9 M. Doi J. Polym. Sci. Polym. Phys. Ed. 1983 21 667. 2o K. Binder J. Chem. Phys. 1983 in press. 21 M. Muthukumar and S. F. Edwards Polymer 1982 23,345. 22 For a review see A. Baumgartner in Monte Carlo Methods in Statistical Physics II,ed. K. Binder (Springer Berlin 1983 in press). 23 A. Baumgartner and K. Binder J. Chem. Phys. 1981 75 2994. 24 M. Bishop D. Ceperley H. L. Frisch and M. H. Kalos J. Chem. Phys. 1982 76,1557. 25 K. Kremer Macromolecules October 1983 in press. 26 A. Baumgartner Proceedings of the Workshop on Dynamics of Macromolecules Santa Barbara 1982 to be published. 27 D. Richter A.Baumgartner K. Binder B. Ewen and J. B. Hayter Phys. Rev. Lett. 1981 47 109; 1982 48 1695. ** B. Ewen B. Stuhn. K. Binder D. Richter and J. B. Hayter to be published. 29 A. Baumgartner J. Chem. Phys.. 1980 72. 871; 1980 73 2489. A. BAUMGARTNER K. KREMER AND K. BINDER 47 30 W. W. Graessley J. Polym. Sci. Polym. Phys. Ed. 1980 18 27. 31 K. Kremer and K. Binder unpublished results. 32 Evidence for reptation in the simulations of lattice chains at higher concentrations is claimed by J. Deutch Phys. Rev. Lett. 1982,49 926. These simulations however are inconclusive as the SAW condition is strictly obeyed for links of different chains while overlaps of a chain onto itself are allowed and thus an artificial mobility of each chain along itself is created such a model leads to reptation trivially but seems rather unphysical. For a more detailed discussion see K. Kremer to be published. 33 J. Klein Nature (London) 1978 271 143; Philos. Mag.,1981 A43 771.