Theoretically Informed Entangled Polymer Simulations : Linear and Non-Linear Rheology of Melts

• Select a chain-end at random, and place the first SL on that chain-end with position Rk(sk). • Create a list of all m segments (a “segment” consists of two successive or bonded beads along the chain) within a distance rcut from the first SL. The cutoff distances rcut is chosen to be larger than bsl . • Select a random position Rnk(snk) (where nk is the partner of the kth SL) along a randomly chosen segment with probability p= exp [ −βusl { Rk(sk)−Rnk(snk) }] /Cnorm. The normalization constant Cnorm is given by Cnorm = ∑l=1 Cl , where Cl is the contribution of the lth segment, explicitly:


Introduction
The general concept behind coarse-grained representations of polymeric liquids is that of universality, the idea that most phenomena at the mesoscale do not depend on microscopic details, but are instead dictated by a few general characteristics of macromolecules, such as their molecular weight and exibility. 1 For example, it is well known that the mean end-to-end distance, R e , exhibits a scaling behavior of the form R e $ N a , where N is the molecular weight or the polymerization index, and where the specic value of the universal exponent a is 1/2 for polymer melts. 1,2A number of universal features also arise in polymer dynamics; for example, the behavior of the viscosity, h, as a function of molecular weight exhibits a scaling behavior of the form [2][3][4] h $ N; N\N c where the value of N c depends on chemical details, but the exponent n is found to have a value of approximately n ¼ 3.4 for a wide range of materials.The quantity N c represents a threshold molecular weight at which the dynamics of polymeric liquids gradually changes from "unentangled" to "entangled".2][3][4] In the so-called tube models, 2,5 originally introduced over four decades ago, the topological constraints on a given chain due to the neighboring chains are replaced by an effective potential, which connes the chain into a tube-like region, thereby coercing the chain to move preferentially along its contour as a result of thermal uctuations.The average time required by a chain to renew its original tube is referred to as the "reptation" time, s rep , and tube models predict a dependence on the molecular weight of the form s rep $ N 3 .As a simple approximation, the viscosity of the polymer can be related to the longest relaxation time, in this case s rep , by h ¼ G * $s rep , where G * is a material-dependent constant.Tube models therefore predict a scaling law of the form h $ N 3 , with a universal exponent that is not too far from that observed in experiments (eqn (1)).
In their simplest form, tube models are known to exhibit a variety of shortcomings; they are, for example, unable to describe the dynamics of polymer blends, the start-up of shear or elongation, or the dynamics of branched polymers, to name a few.Over the last decade, important advances have been made to address such shortcomings.Recent state-of-the-art descriptions are not limited to reptation as the only relaxation mechanism, but consider additional modes, including (a) contour length uctuations (CLFs) and (b) constraint release (CR).Contour length uctuations are associated with changes of the tube length due to accordion-like motions of the chain inside a tube. 6Constraint release is a multi-chain effect; the molecules surrounding the probe chain also undergo relaxations, thereby leading to reorganization of the probe's tube.CLF and CR are critically important in polydisperse samples, where, for example, the constraints imposed on longer chains by surrounding shorter chains have short life times, and facilitate the relaxation of the overall blend in a non-trivial manner.The addition of CLF and CR to the tube model by various authors has enabled prediction of linear viscoelastic properties for different classes of homopolymer systems, and available implementations have been discussed in a recent review. 7 particularly useful variant of tube models for the rheology of polymeric liquids relies on the concept of slip links.Slip links, originally introduced by Edwards et al. 8,9 to represent the effective constraints on a chain due to entanglements, tether the chain contour to its primitive path, thereby constraining its motion.Different implementations of slip-link based models by Schieber et al. 10,11 and Likhtman et al. 12,13 have established the general usefulness of the slip-link concept for a wide variety of systems.More recently, Müller and Daoulas 14 have taken an additional step and combined the single-chain-in-mean-eld (SCMF) approach 15 with slip links, whose temporal evolution was implemented using a Monte Carlo approach.In all such approaches, an individual chain is tethered to the slip links, and their extension to complex ows or inhomogeneous systems is unclear.
7][18][19] These primitive paths interact in a pairwise manner via slip-links, whose temporal evolution is governed by a force balance on the entanglement points.The monomers can slide through the sliplinks at a rate governed by the distribution of local tension in the chain.Such models have also been shown to provide good agreement with experimental linear response data for homopolymer melts.In such approaches, however, the molecular structure of polymeric molecules is not resolved, and extension to nanostructured or inhomogeneous systems is challenging. 40urthermore, the nodes (entanglements) agglomerate into clusters; one must therefore introduce arbitrary repulsive interactions between them, without which the model is unable to capture the equilibrium properties of homopolymers. 20ast implementations of slip-link models of primitive-path network models have generally been limited to melts of pure homopolymers.Notable exceptions are provided by the studies of Schieber and co-workers and Masubuchi and co-workers, who also examined the linear response of various bidisperse homopolymer blends and found good agreement with experiments.Available models have not been extended to nanostructured, microphase separated polymeric materials.In this work we present an approach that builds upon earlier versions of slip-link representations of dynamics and that incorporates the thermodynamics of the system, thereby enabling description of linear and non-linear dynamics in materials having complex molecular architectures and undergoing phase transitions or microphase separation.The validity of the proposed model is demonstrated by examining the linear and non-linear response of pure homopolymers and bidisperse blends for which experimental data are available, and by examining the anisotropic diffusion that arises in microphase separated diblock copolymers.4][25] The slip-links are introduced at the two-chain level, as uctuating interactions between pairs of polymeric molecules, thereby enabling a facile description of dynamics in inhomogeneous systems.It is found that for the systems considered in this work, the proposed approach provides a good description of experimental observations.
It is important to note that the purpose of this work is not to provide an exhaustive description of available rheological data for entangled polymers, but to demonstrate instead the ability of a coarse-grained description typically used in eld theory to describe the dynamics of homogeneous and inhomogeneous polymeric materials.Our results indicate that the proposed model can do that, and suggest that additional studies with a variety of microstructured systems may be warranted.

Model and methods
Without loss of generality, we present the method in the context of an AB diblock copolymer.The melt consists of n chain molecules in a volume V at temperature T. The macromolecules are represented by exible linear chains described by the discretized Gaussian chain model, in which each chain is composed of N beads.The position of the s th bead in the i th chain is given by r i (s).Thus, the intra-molecular interactions acting along the chain's backbone are given by bH where N ¼ (r o R e 3 /N) 2 , where r o is the bead number density and R e 2 ¼ (N À 1)b 2 is the mean squared end-to-end distance of an ideal chain.The repulsion between unlike monomers, represented by the rst term in the functional, is characterized by the Flory-Huggins parameter c.The second term, given by Helfand's quadratic approximation, 27,28 ascribes a nite compressibility to the melt, which is related directly to kN. 15,29 Note that ffiffiffiffi N p , the interdigitation number, represents the average number of molecules that a chain interacts with, and is proportional to the molecular weight of the polymer.9][30][31] Intermolecular interactions are then expressed as pair-wise interactions, bH nb ¼ P i\ j U ij , where the sum runs over all particles in the melt. 28,29With this formalism, the intra-and inter-molecular forces acting on each particle can be obtained directly, f b and f nb respectively.As expected, the resulting inter-molecular interactions are so i.e., the interaction energy remains nite in the limit r / 0, where r is the distance between particles.Therefore, entangled dynamics cannot be reproduced by this model in its bare form.As alluded to earlier, the effects of entanglements are introduced through the use of slip links (SLs), as illustrated schematically in Fig. 1.For conciseness, we refer to the proposed implementation as "Theoretically Informed Entangled Polymer Simulations" (TIEPOS).In that implementation, we assume that an entanglement is a binary interaction chains.Thus the SLs evolve in pairs, consisting of two rings that slip along different chain contours, and are connected by a spring.The positions of the rings fR ik g Nsl k¼1 , where N sl is the total number of slip links in the system, depends parametrically on a continuous scalar variable s k ˛[0, N À 1] on polymer i k , i.e.R i k (s k ).Note that slip-links can "live" on different chains.The rings move in straight lines between beads.To include these topological interactions at the same level of description as that of the polymer chains in the model, we use a harmonic interaction between the rings: bu sl ðs k ; where n k is the partner of ring k and b sl characterizes the strength of the spring.N sl and b sl are parameters to be determined by comparison to experimental data.Given the characteristics of a specic polymer (molecular weight, c parameter, etc.) and a specic contour discretization N, N sl can be inferred directly from knowledge of the molecular weight between entanglements for the particular polymer melt of interest.Thus, if N o sl is the average number of entanglements in a chain, N sl ¼ nN o sl .For the strength of the spring we choose b sl ¼ b.
The temporal evolution of the material is governed by the following stochastic differential equations for beads, j ¼ 1,., nN, and for the rings, k ¼ 1,., N sl .In eqn (3), f sl represents the forces acting on the beads due to the slip links, g is the associated friction coefficient, and x represents the stochastic inuence of the coarse-grained microscopic degrees of freedom and satises the uctuation-dissipation relationship: , where a, b ¼ x, y, z.In eqn (4), h sl k is the projection along the respective bond of the slip-link force and g sl is the friction coefficient of the rings.Parameter z is a random force with zero mean that satises the uctuationdissipation theorem.Note that we have chosen to implement a simple Brownian dynamics approach, but other approaches, such as DPD, MD, etc., could be equally implemented.
SLs restrict the lateral motion of the chain with respect to the axis of the backbone, and favor the movement along the chain contour.The chain contour length naturally uctuates due to the Gaussian chain representation adopted here.Thus, two relaxation mechanisms are included in the model, namely reptation and CLF.The other main relaxation mechanism, constraint release, is implemented as follows: rst, to initialize a new SL pair, a chain-end is randomly chosen and the rst ring of the SL is attached to that position.The partner ring is then attached to a nearby chain in such a way as to generate a Boltzmann distribution of links (see ESI †).Once initialized, the dynamics makes it possible for one ring to leave its chain, i.e. s k ; [0, N À 1]; when this occurs, the ring and its corresponding partner could be destroyed (see ESI †); when that happens, a new SL pair is created as described above.This process mimics the disentanglement at the chain-ends and the constraint release.
The TIEPOS approach proposed here should be contrasted with primitive chain network models, where the nodes are the entities that move in real space, and polymer chains are displaced because the rings pull them over.In the proposed TIE-POS approach, the chains evolve on their own under the constraints provided by the SL.Also note that SLs introduce an additional effective attractive potential between the chains.The functional form of that potential can be shown to be of the form u eff ðrÞ ¼ Àzk B Te ÀuslðrÞ=kBT , where z is the activity of the slip links and r denotes the distance between beads. 32In melts, which are the focus of this work, the structure and thermodynamic properties of the polymer are found to be largely unaltered by the inclusion of the slip links.The thermodynamic behavior is dominated by the incompatibility (cN) and compressibility (kN) contributions to the energy (eqn (2)); the compressibility term is high enough to avoid the clustering presented in primitive chain network models.Note that in dilute solutions one must include the potential v eff ¼ Àu eff in order to cancel the effect of the slip links on thermodynamics.Dilute solutions, however, are unlikely to be entangled and, for such systems, the formalism proposed here is of limited relevance.Unless otherwise noted, all of the calculations reported here were performed without v eff .In addition to verifying that the thermodynamics are unaltered by the inclusion of slip-links, we have also veried that the distribution of links along the contour of the chain is uniform.The corresponding results are shown in Fig. S1 of the ESI.† For the particular case of microphase separated copolymers, it remains to be veried whether the slip-links omitting v eff alter drastically the equilibrium phase behavior or not.Our simulations suggest that the SL effect is negligible, but a more extensive analysis will be presented in a future publication.

Results
In the absence of SLs, the model described above reproduces Rouse dynamics (data not shown).The critical issue addressed here is whether the simple implementation of slip links described above reproduces the behavior of entangled polymers.For concreteness, in our simulations we have chosen parameters that are representative of polystyrene (PS) melts, a commonly used polymer for which extensive experimental data are available.We consider two different monodisperse samples: (i) PS with a molecular weight w ¼ 103 kg mol À1 (PS100) and (ii) PS with a molecular weight M w ¼ 176.7 kg mol À1 (PS177).Experimental loss and storage modulus data for these materials are available in the literature from Nielsen et al. 33 and Maier et al. 34 The stochastic differential equations (eqn ( 3) and ( 4)) are solved using the Euler algorithm.The time is rescaled by the characteristic diffusion time of a system with only random forces, s 0 ¼ l 0 2 /D 0 (s 0 is the diffusion time of a bead and l 0 is a reference length scale chosen as l 0 ¼ R e ¼ 1), where D 0 ¼ k B T/g is the diffusion constant of unconnected monomeric units under these conditions; the rescaled time is given by t * ¼ t/s 0 .In these units, the time step used to integrate the equations of motion is given by Dt * ¼ 10 À4 ; the energy is measured in units of k B T. The bead number density is xed at a reference value of r o ¼ 4096/R eo 3 , where R eo is the reference scale length of a PS homopolymer with M w ¼ 103 kg mol À1 represented by 32 beads.
The parameters used to generate our results were c ¼ 0, k ¼ 1.5625, and g sl /g ¼ 0.1.The molecular weight between entanglements used here was M e ¼ 15.65 kg mol À1 . 41Table 1 summarizes the parameters used in our calculations; the same set of parameters was used for the monodisperse systems and for the blends (see below).We compute the loss and storage moduli as indicated in the ESI.† The comparison to experimental data is accomplished by scaling the frequency u and the modulus G 0 ; the two tting parameters in the model are determined from the experimental G 0 and G 00 curves at the point where they intersect each other.The scale factor obtained in this manner is then applied to all numerical data.Results for the two different monodisperse samples considered here are shown in Fig. 2 and 3.The number of entanglements per chain was N o sl ¼ 7 for the PS100 and N o sl ¼ 11 for PS177.It can be seen that the results from simulations are in quantitative agreement with experimental measurements over almost ve decades in the frequency domain.Discrepancies start to appear at high frequencies, which correspond to fast processes that cannot be captured by a coarse grained model (which omits truly microscopic, atomistic processes).Additional results for the scaling of diffusion and the longest relaxation time with molecular weight are provided in the ESI (Fig. S2 †).There is good agreement between our predictions, experiment and theory.
As mentioned above, a critical test of the model is whether it is able to predict the behavior of more complex materials, such as bidisperse samples comprised of short and long molecular weight chains.For highly asymmetric components, the CLF and CR relaxation mechanisms become considerably more important than for monodisperse samples.For concreteness, we consider bidisperse PS samples with molecular weight M w ¼ 36.3 kg mol À1 (PS36) and M w ¼ 407 kg mol À1 (PS407), which have been characterized experimentally. 35We consider two different weight fractions of the low molecular weight component, namely w low ¼ 0.8 and w low ¼ 0.6.For these systems, the entanglements per chain were N o sl ¼ 26 for the large chains (PS407) and N o sl ¼ 2 for the short chains (PS36).Results for the storage and loss moduli are shown in Fig. 4. The agreement between simulations and experimental data is quantitative, suggesting that the main relaxation mechanisms that arise in entangled polymers are correctly captured by the proposed TIEPOS approach.Fig. S2 of the ESI † also shows results for homopolymers that indicate that the viscosity scales with the 3.4 exponent of molecular weight.
One of the key features of the TIEPOS formalism adopted here is that the Hamiltonian can be used to describe  View Article Online nanostructured macromolecular systems, such as microphase separated block polymers.We apply it here to the equilibrium dynamical behavior of a symmetric diblock copolymer below the order-disorder transition.For the calculations employed here we use kN ¼ 35, cN ¼ 25 and ffiffiffiffi N p ¼ 110.Such parameters have been shown to provide a quantitative description of PS-PMMA diblock copolymers with molecular weight 104 kg mol À1 .In units of the end-to-end distance, the lamellar period in the bulk corresponds to L 0 ¼ 1.66R e .The chains are discretized into N ¼ 16 + 16 beads.We use periodic boundary conditions in the three spatial directions.For completeness, we compare the dynamical behavior of this nanostructured polymeric material in the limit of unentangled (Rouse) chains, and when the TIEPOS approach is used to describe entanglements.We describe the entangled state by incorporating 6 SLs per chain.For simplicity, we assume that the distribution of entanglements in a lamellar morphology is similar to that of a disordered melt; we recognize that entanglements could exhibit different distributions for different morphologies but, given the limited information that is currently available for entangled copolymer dynamics, we feel that such an assumption is warranted.
Due to the lamellar structure of the material, it is expected that the behavior in the two orthogonal directions, parallel and perpendicular to the interfaces, will be different.We therefore compute the mean square displacements for the beads g 1 (t) and the centers of mass of the chains g 3 (t) along the parallel and perpendicular directions.Barrat and Fredrickson have predicted that, for unentangled (Rouse) chains, the motions parallel and perpendicular to the interfaces are decoupled. 36ig. 5 shows results for Rouse chains; it can be observed that the parallel motion is not affected by the perpendicular dynamics, thereby conrming Barrat and Fredrickson's predictions.The perpendicular dynamics exhibits a plateau, which corresponds to the enthalpic penalty associated with chains "crossing" the lamellar domains.Such crossing events are slow, and decrease the dynamics.By computing the diffusion constants associated with motion in the parallel and perpendicular directions we arrive at a ratio D k D t z80 for unentangled chains.Barrat and Fredrickson 36 predicted a diffusion coefficient ratio of O (10), with an exponential dependence on cN, again consistent with our own ndings.Fig. 5 shows the results when entanglements are introduced into the model.In this case the perpendicular and parallel motions are no longer decoupled. 14The motion of the center of mass as a function of time in the parallel direction is no longer linear (on a log-log scale), as it was for unentangled chains (see Fig. 5a), but instead it exhibits a weak plateau, which coincides with that observed for motion in the perpendicular direction.The asymmetry in the perpendicular and parallel diffusion coefficients is much less pronounced than in the unentangled case.The parallel motion continues to be faster than the perpendicular, but the diffusion coefficient ratio is now well below that observed for unentangled chains.In a set of seminal experiments, Lodge and Dalvi measured this ratio for moderately entangled diblock copolymers. 37,38They found a value of approximately 4; which is again consistent with our entangled simulations, and serves to emphasize the plausibility of the model and approach proposed here.Note that in this rst implementation of the formalism we have not attempted to provide a quantitative comparison with experiments; in addition to assuming that the entanglements follow the same distribution as in the disordered homopolymer melt, we have also assumed that the friction coefficients of different species in different domains are the same.These assumptions are likely to require revisions once we examine specic systems.The main point here, however, is that the formalism proposed here will enable systematic investigations of such assumptions.It is of interest to compare the TIEPOS method to that proposed by Masubuchi et al., 21,22 who also considered the dynamics of an entangled diblock copolymer.In their work, these proposed a free energy that is an explicit function of the entanglements.In our approach, the underlying free energy functional depends on local composition, and it is the same as that employed in eld theoretic treatments of block polymers.At this point, we have demonstrated that the proposed approach correctly describes the linear dynamics of melts.We conclude this study by also examining its ability to describe the non-linear rheology of polymer melts under shear ow.We use Lees-Edwards boundary conditions along the gradient direction (y-direction); the ow proceeds along the x-direction.A ow eld is imposed in the differential equations that reproduces the linear velocity prole in shear ows.The parameters are chosen such that the model represents a monodisperse PS melt with a molecular weight M w ¼ 100 kg mol À1 .We compute the steady shear viscosity h ¼ s xy / _ g for different shear rates _ g.The results are shown in Fig. 6.The rst point to note is that the model reproduces the shear thinning behavior found experimentally.Second, viscosity follows a power law scaling as h $ _ g À0.75 in the shear-thinning regime.This exponent is similar to that observed experimentally (around À0.71). 39Predictions for elongation (data not shown) are also consistent with experiments.Additional results for model polymers of different molecular weights are included in the ESI (Fig. S3 †); the scaling with molecular weight predicted by the model is consistent with experimental data and theory.

Conclusions
In summary, by combining the concept of slip-links with a widely used theoretically informed coarse-grained polymer model, we have arrived at a simple, yet useful formalism for prediction of the linear and non-linear rheology of polymeric melts.The validity of the model has been established by examining its predictions for the storage and loss moduli of polystyrene blends, by analyzing the diffusion of polymer molecules in a microphase separated symmetric block copolymer, and by determining the viscosity of a homopolymer melt in the shear thinning regime.In all cases, the agreement with experiments appears to be quantitative, serving to validate the overall structure of the proposed approach and suggesting that it may have general applicability to more complex materials, including nanocomposites, concentrated solutions, and molecules having a more complex molecular architecture.
where b 2 is the mean squared bond length of an ideal chain, b À1 ¼ k B T, and k B is the Boltzmann constant.Inter-molecular interactions are represented as a function of the local densities f A,B (r):

Fig. 1
Fig. 1 Schematic representation of the introduction of slip links in a many-chain model.

Fig. 5 Fig. 6
Fig.5Temporal evolution of the mean square displacements for beads g 1 (t) and centers of mass of the chains, g 3 (t), computed along the parallel and perpendicular directions to the lamellar interfaces: (left) unentangled chains and (right) entangled chains.

Table 1
Parameters characterizing the homopolymer melts