Article pubs.acs.org/Macromolecules
Finding Entanglement Points in Simulated Polymer Melts Jing Cao, Jian Qin, and Scott T. Milner* Department of Chemical Engineering, The Pennsylvania State University, University Park, Pennsylvania 16802, United States ABSTRACT: There are two complementary views of entanglement in polymer melts: a continuous tube confining transverse motions of a given chain and discrete binary entanglement points with neighboring chains. The discrete view is implicit in chain-shrinking studies of simulated entangled melts, while our recent isoconfigurational ensemble (ICE) technique determines a continuous primitive path and confining potential. Here, we reconcile these competing views by identifying points of close contact along the ICE primitive path in an entangled ring polymer melt. We set a minimum distance criterion such that the number of close contacts is equal to the number of entanglement strands, which gives on average about two nearby entanglement points acting on any primitive path segment, a sensible result. We show that both the continuous primitive path and the discrete close contacts are long-lived for ring melts, for which topological confinement effects are expected to be permanent.
■
INTRODUCTION
In a second, complementary view, the tube is represented as a continuous path in space, called the primitive path, to which the given chain is bound by a confining potential. The primitive path is assumed to be a semiflexible random walk with persistence length of order a and total contour length of (N/ Ne)a. Correspondingly, the confining potential is assumed to confine transverse fluctuations of the chain to distances of order a from the primitive path. The tube diameter a and the corresponding value of Ne are material parameters, which depend on monomer architecture and solution concentration, but are independent of details of chain architecture such as chain length and the presence of long branches. There exists both a well-established correlation for determining a from chain dimensions (persistence length and monomer volume)12 and a theoretical understanding of the physical origin of this correlation.13 The discrete view of entanglements can be regarded as a coarse-grained version of the continuous tube, in which the smooth primitive path is replaced by a sequence of straight steps of length a between pointlike constraints. This complementarity is sometimes computationally convenient2 but invites the questions: Which view is more microscopically realistic? Is the tube really continuous? Do entanglement points really exist? Simulations are well-suited to answering these questions, since within a simulation we can visualize the conformations and trajectories of individual chains and their neighbors. So how should we go about finding the tube, or entanglement points, in simulations? One approach that has been fruitfully employed to investigate entanglements in simulated polymer melts is the chain-shrinking method, developed by Everaers et al.14 and
Contemporary theories for the flow behavior of entangled melts and solutions are based on the idea of a tube confining the transverse motion of a given chain, as a consequence of the uncrossability of neighboring chains.1,2 Starting from the assumption of a confining tube, theories of stress relaxation after a step strain focus on mechanisms by which a chain escapes its tube and explores new conformations, forgetting the anisotropic conformations induced by the step strain. For linear chains, these mechanisms include contour length fluctuations, constraint release, and reptation. For branched chains (stars, Hpolymers, and the like) chain segments with free ends explore new tubes by arm retraction.2−7 Employed together, these mechanisms give a good quantitative account of stress relaxation in a wide variety of entangled melts and solutions.7−11 The tube itself is an intuitively appealing concept that remains somewhat elusive, both in experiments and in simulations. There are two complementary views of the confining effects of entanglement of a given chain by neighboring chains, both of which have been fruitfully employed in theoretical calculations. In one view, the tube is regarded as a sequence of entanglement points of a given chain with neighboring chains. These entanglement points are sometimes modeled as “sliplinks”, which themselves may be either fixed in space or attached by Gaussian springs to nearby fixed points. The number of monomers along a chain between successive entanglement points is Ne, the entanglement molecular weight. The mean-square end-to-end distance of such an entanglement strand is the tube diameter a. The trajectory of a given chain in an equilibrated melt is taken as a discrete random walk of steps between successive entanglement points, with N/Ne steps of step length a. © XXXX American Chemical Society
Received: May 17, 2014 Revised: December 1, 2014
A
dx.doi.org/10.1021/ma5010315 | Macromolecules XXXX, XXX, XXX−XXX
Macromolecules
Article
Larson et al.15 In the chain-shrinking method, the primitive paths of simulated entangled chains are obtained by fixing the chain ends in space and switching off intramolecular interactions and then minimizing the energy. In this way, system excluded volume is greatly reduced, but chains still cannot cross (because interchain interactions are preserved). The resulting chain conformations are piecewise linear paths between binary entanglement points, at which two chains cannot cross. Similar chain-shrinking methods were developed by Kröger16 and Theodorou17 using a geometric approach to straighten chain contours without allowing chains to cross. Using chain-shrinking methods, the mean number of monomers between entanglement points can be measured, from which Ne and hence a can be determined. Of course, the chain-shrinking process destroys the local structure of the melt, so that details of the tube and its confining potential cannot be observed. Indeed, one may say that the chain-shrinking approach subtly reinforces a discrete view of binary entanglements between chain strands, without actually giving support to the existence of binary entanglements in actual melts, since chain shrinking inevitably produces a configuration of straight paths between crossing points. A different approach to visualizing the tube is provided by the recently developed isoconfigurational ensemble (ICE) method, which gives a noninvasive view of chain fluctuations within the tube.18−20 In the ICE method, multiple short molecular dynamics simulations are performed starting with the same entangled configuration, with different initial velocities. The primitive path is obtained by averaging the chain conformations over the different trajectories and over a short averaging time τa, typically of order the entanglement strand relaxation time τe. The tube confining potential can be obtained from the shape of the “cloud” of monomer positions surrounding the primitive path. The primitive path so obtained is smooth, with a persistence length of order the tube diameter; the confining potential is found to be harmonic, with a width of order a. The values of a and hence Ne derived in this way are consistent with those obtained from chain-shrinking methods as well as values obtained from a completely different approach of measuring the topological entropy of the ring melt.21 Just as the chain-shrinking methods confer a subtle bias toward regarding entanglement points as discrete, so the ICE method implicitly assumes a view of the tube as continuous. The ICE method is mean field in the sense that an average is taken over all the ways the chain could fluctuate, resulting in a smooth primitive path that is nearly local and nearly instantaneous, smoothed on the length scale of the tube diameter a and time scale of the entanglement strand relaxation time τe. In this work, we will present a simple method to identify entanglement points “nondestructively”, without resorting to the invasive methods of chain shrinking, by locating points of close approach of the ICE primitive path. In a recent paper, Likhtman has presented a more sophisticated algorithm for defining the primitive path, called the “tube axis algorithm”, which he uses to investigate the properties of entanglements.20 Here we briefly summarize his approach and main conclusions. Reference 20 defines a time-averaged primitive path R(s) from a single long time trajectory of entangled chains (no isoconfigurational averaging is performed), in which many chain configuration “snapshots” ri(s) (i = 1, ..., M) are collected. The primitive path R(s) is defined as
R(s) = 1/M ∑ ri(si(s)) i
(1)
in which si(s) is a monotonic mapping of the arc length coordinate of the ith configuration to the primitive path. The mappings si(s) are chosen to minimize the sum of the squared distances between the remapped configurations ri(si(s)) and the primitive path R(s). The minimization is performed by means of an equivalent slip-spring model, for which simulated annealing is used to carry out the minimization. The result of this procedure is a smooth primitive path R(s), which as a result of the remappings is insensitive to motion of the chain within the tube. Likhtman applies the tube axis algorithm to analyze simulations of melts of multiple entangled rings, introduced by our group to study the tube and its properties (see below), using the standard Kremer−Grest model. He presents results for three systems: (a) 100 rings of N = 256 beads, (b) 150 rings of N = 512 beads, and (c) 240 rings of N = 1024 beads. All three systems were run for a total time of (1−3) × 106τLJ, where τLJ is the Lennard-Jones time (see below). In ref 20, the entire time trajectory is used to define the primitive path; thus, the path is time-averaged over a very long time, many times longer than the curvilinear Rouse time (stretch relaxation time) of the ring in its tube. Thus, the primitive paths in that work are not “nearly instantaneous”, as are the ICE primitive paths we use here, which are averaged only over τe or so. One fundamental quantity to emerge from the tube axis algorithm is the tube radius, defined as the average distance of monomers from the primitive path. Reference 20 reports a tube radius of about 2.5σ, quite consistent with our previous results using ICE averaging.19 Reference 20 emphasizes the curvature of the primitive path as a sensitive indicator of interactions with other chains. However, entanglement points are most simply identified in terms of points of close approach between different parts of the primitive paths (an approach we also use in this work). Entanglement points on a given ring are found to exhibit repulsive correlations, with a range of about ten monomers. A heuristic value of close approach threshold gives Ne ≈ 35−40 as the average distance between contacts, a value substantially smaller than our best estimates of Ne by topological methods and ICE analysis.19,21,22 A certain number of configurations are reported in ref 20 in which three or more primitive path strands coentangle, consistent with our previous investigations of topological changes resulting from chain crossing in entangled melts.23 However, the fraction of these is estimated to be in the range of 20−30%, whereas we observed such configurations to be rather rare (a few percent) in our chain crossing studies. Reference 20 describes the purpose of the work reported there as “to understand what entanglements and tube axes are, rather than to count them”. In contrast, our emphasis in the present work is to identify entanglement points along the primitive path by a simple heuristic and to answer two basic questions: (1) Are there the right number of entanglement points? We expect to have of order one entanglement point per entanglement strand. (2) Do the entanglement points identified by the heuristic persist in time? A melt of entangled rings is ideal for answering both of these questions. Indeed, rings are particularly well suited to simulation studies of the tube. For entangled rings, the tube B
dx.doi.org/10.1021/ma5010315 | Macromolecules XXXX, XXX, XXX−XXX
Macromolecules
Article
Monte Carlo chain crossing moves used for topological equilibration in the section Isoconfiguration Ensemble (ICE) Primitive Path. There, we briefly describe our method for noninvasively determining the primitive path. We show how our primitive path depends on the averaging time τa and select a convenient averaging time for our present study. In the section Close Contacts on the Primitive Path, we present our approach for identifying persistent entanglement points, which is based on the evaluation of the squared distance r2(i,j) between points i and j on the primitive path. The landscape of this function in the i−j plane has many minima, some of which correspond to close approaches between primitive path segments. We choose a cutoff distance rc to select from among these minima those corresponding to entanglement points, with the requirement that we obtain on average N/Ne points. Having identified the entanglement points by this heuristic but well-motivated procedure, we can ask whether the entanglement points persist in time. In the section Persistence of Entanglement Points, we investigate this question qualitatively by imaging the world lines of entanglement points25 as well as quantitatively by computing their timedependent correlation function. By both means, we find that indeed the entanglement points we identify do persist, as we expect them to do, since the entanglement state of a ring melt is permanent. We do observe some “chattering” on and off of entanglement points, as the ring melt undergoes various thermal fluctuations consistent with the given topological state.
is expected to be permanent. With no free ends to relax, processes such as tube relaxation due to contour length fluctuations, constraint release, and reptation are completely suppressed. We have studied melts of long entangled rings entangled with their periodic images as a proxy for the local state of entanglement of a melt of long entangled chains. For many purposes we prefer to study a single long self-entangled ring, one reason being that finite-length corrections are smaller for a given system size. For such simulations, we must equilibrate the ring or rings over the ensemble of topological states, which are permanent unless we introduce some way for the chain to cross itself. To topologically equilibrate ring melts, we perform hybrid MD/ Monte Carlo simulations, in which we introduce various chaincrossing moves.21 These lead to what amounts to constraintrelease Rouse motion of the rings, as the chain segments infrequently cross. We monitor the Rouse mode fluctuations of the rings facilitated by the chain crossings to verify that we have equilibrated the ring conformations with respect to entanglement. Then, to study the effects of chain uncrossability on the local motion of the chains, we perform various MD simulations with the chain-crossing moves turned off. We have used this approach to prepare topologically equilibrated self-entangled ring melts, which we have then studied with ICE methods to obtain tube properties.19,21,24 Our strategy for identifying entanglement points is based on the ICE primitive path. In brief, we look for points of close approach between two distinct segments along the path. This effectively introduces some preaveraging into the search for entanglement points, since the ICE primitive path is averaged over multiple trajectories starting from the same configuration, as well as time-averaged over a short averaging time τa. The preaveraging is useful because it leads to a smooth primitive path that likewise varies smoothly in time. The ICE primitive path varies with the isoconfigurational averaging time τa. As τa increases, the primitive path becomes increasingly smooth, as more slowly relaxing, longer-wavelength contour-length fluctuations within the tube contribute to the time average. An artifact of this smoothing, called “corner cutting”, causes two entangled portions of primitive path to appear to cross through each other as the averaging time is increased. The tube axis algorithm of ref 20 avoids this artifact, as a result of the remapping applied to the arclength coordinates in defining the primitive path. Here, we take advantage of the corner cutting artifact, choosing the averaging time τa such that tightly entangled portions of the primitive path appear to nearly intersect. The convenient choice of τa turns out to be near the entanglement strand relaxation time τe. This is a pleasing result in the sense that the physical size of a tight entanglement point might be expected a priori to be about one tube diameter. Having determined the ICE primitive path, we identify entanglement points by selecting a cutoff distance rc for a close approach between distinct portions of the primitive path. We select rc heuristically, so that on average we find about N/Ne entanglement points. It turns out that the rc value that gives about N/Ne entanglement points also gives the pleasing result that any given primitive path strand participates in or is influenced by about two contact points on average (a statement we will make more precise in the section Close Contacts on the Primitive Path). The balance of this paper is organized as follows. We summarize our molecular dynamics simulation model and the
■
ISOCONFIGURATIONAL ENSEMBLE (ICE) PRIMITIVE PATH For our simulations, we use a coarse-grained bead−spring model for polymer rings.26 The Hamiltonian includes a purely repulsive Lennard-Jones pair potential between nonbonded monomers Up(r ) = 4ε((σ /r )12 − (σ /r )6 + 1/4)
(2)
with ε = kBT, and a harmonic potential between bonded neighbors Ub(r ) = (1/2)κ(r /σ − 1)2
(3)
corresponding to a finite length spring (rest length σ) with a spring constant κ = 400kBT. The nonbonded potential is smoothly cut off at short distances for interaction energies above 80kBT to avoid excessively large forces. The bead number density ρ is set to be 0.7σ−3, to represent a melt. The average bond length and the statistical segment length b are found to be σ and 1.414σ, respectively.26 We simulate a single long self-entangled ring of N = 800 monomer beads in 3d periodic boundary conditions, as a proxy for an entangled melt of long linear chains. The corresponding linear dimension L of the system is L = 10.46σ. The initial configuration of the ring is created by generating a random walk of N steps, computing its end-to-end vector R, and then adding −R/N to all N bonds so that the ring closes. The configuration is prepared for simulation by (1) using random bead displacement MC move to eliminate bead overlaps and (2) using a combination of hybrid MD/MC and rebridging moves to equilibrate the pressure and topology. In a hybrid MD/MC move, a short MD run is used to generate a trial configuration, to which the Metropolis acceptance criterion is applied. C
dx.doi.org/10.1021/ma5010315 | Macromolecules XXXX, XXX, XXX−XXX
Macromolecules
Article
We topologically equilibrate our self-entangled ring by introducing a variety of Monte Carlo moves that permit the ring to locally cross itself.21 Without these moves, the topological state of the ring in our simulations is permanent. We monitor the equilibration by observing the resulting constraint-release Rouse motion of the ring, which permits us to explore the full range of Rouse fluctuations without regard to entanglement.23 The fundamental microscopic time scale of our bead−spring MD simulations is the Lennard-Jones time τLJ, defined by
τLJ = (mσ 2/ε)1/2
(4)
in which m is the bead mass and σ and ε are the Lennard-Jones parameters. (Particles with kinetic energy ε travel a distance of order σ in the time τLJ; hence, τLJ is the time scale for collisions to occur in a moderately dense Lennard-Jones fluid, for which ε is of order kT.) In this work as for entangled polymer rheology generally, the fundamental “mesoscopic” time scale is τe, the Rouse time of an entanglement strand. In a forthcoming publication, we have determined τe for our system of a self-entangled ring polymer melt by comparing simulation results for the monomer meansquare displacement versus time to an analytical theory.24 Roughly speaking, the diffusive motion of monomers ⟨Δr2(t)⟩ crosses over from unhindered three-dimensional Rouse motion to one-dimensional, curvilinear Rouse motion along the primitive path at the time scale τe and length scale a. From the detailed analysis of ⟨Δr2(t)⟩ we conclude that τe ≈ 2000τLJ for our simulation model. To obtain the primitive path of our simulated chains, we use isoconfigurational ensemble (ICE) averaging. Starting from the same initial configuration, a series of independent MD simulations are performed (with chain-crossing moves turned off, so the chains cannot cross), with different sets of random thermalized initial velocities for the monomers. In our work, typically 100 independent simulations are generated. In each simulation, time trajectories are recorded, with an interval between “snapshots” of 50τLJ. The primitive path is obtained by averaging the positions of each monomer in the resulting trajectories over the isoconfigurational ensemble (i.e., over the choice of initial velocities) and over a short averaging time τa. The averaging time τa is typically taken to be of order the entanglement strand Rouse time τe. On this time scale, a chain segment first explores the confining tube, which arises from entanglement with nearby chains. The ICE procedure results in a smooth primitive path that is as local as possible in space and time (the primitive path cannot be defined on time scales much shorter than τe and hence on length scales much shorter than a). Figure 1 shows the primitive path resulting from ICE averaging of one self-entangled ring of N = 800 monomers, for different averaging times. The chain configuration and primitive paths have been “unfolded” from the original periodic simulation volume (cubic box in the center of the figure), into smooth closed trajectories. Evidently, the unfolded ring is much larger than the simulation box, and the ring entangles with its periodic images. The ICE primitive paths for different averaging times τa are shown (green: 0.2τe; blue: 0.8τe; purple: 2τe). The ICE primitive path depends to some extent on the choice of averaging time. With increasing τa, the primitive path is smoother, with shorter contour length, as fluctuation modes
Figure 1. Initial chain conformation (red path and points) and primitive paths, “unfolded” from the periodic box (black outline), for different averaging times (green: τa = 0.2τe; blue: τa = 0.8τe; purple: τa = 2τe).
with relaxation times shorter than τa are averaged over. Because of the simple form of averaging used in our method, in which the primitive path position of a given monomer is taken to be its ICE average position, an artifact called “corner cutting” appears in the primitive path conformations as the averaging time is increased. Consider a monomer at some point along the primitive path at which the path curves sharply around some other entangled strand. This monomer fluctuates back and forth inside the tube, resulting in a cloud of positions that becomes increasingly arcshaped as τa is increased. Correspondingly, the average position of the cloud migrates toward the center of curvature of the true primitive path. The result is that a pair of tightly entangled primitive path segments appear to pass through each other, as the averaging time increases. This behavior is evident in Figure 2, in which the entire primitive path is displayed within the periodic system volume, for different averaging times. Two entangled primitive path strands are highlighted (red and green). For short averaging time τa = 0.1τe (top left), the two strands are clearly linked and curved away from the entanglement point in response to tension forces along the primitive path. As the averaging time increases from τa = 0.2τe (top right) to τa = 0.4τe (bottom left), the two strands of primitive path appear to pass through each other. The effect is more pronounced for τa = 0.8τe (bottom right) and longer averaging times. Our main concern in this paper is to develop a way to identify discrete entanglement points along the primitive path. In the next section, we will exploit the tendency of mutually entangled primitive path strands to pass through each other as the averaging time increases, to identify entanglement points. In brief, we will choose τa to be around 0.2τe, precisely because for averaging times in this range, tightly entangled primitive path strands appear to very nearly intersect each other. We then look D
dx.doi.org/10.1021/ma5010315 | Macromolecules XXXX, XXX, XXX−XXX
Macromolecules
Article
entanglement point. To find entanglement points based on the raw, unprocessed chain trajectories is very challenging because the location of an entanglement “point” necessarily must be uncertain to within a or so in space and to within Ne or so along the entangled chain strands. Thus, many monomers somehow participate in an entanglement point, and it is somehow unnatural to try to localize the entanglement point to a single pair of monomers. Our ICE primitive path provides a preaveraging of the complete chain trajectories which is well suited to identifying entanglement points. As we have seen, the resulting primitive paths are smooth on the scale of the tube diameter. Furthermore, we can choose the averaging time τa so that tightly entangled primitive path strands appear to nearly intersect each other. This motivates us to define entanglement points in terms of the map of spatial distances r(i,j) between points i and j on the primitive path. The distance map r(i,j) for a given primitive path R(i) is simply the 2d square array of distances |R(i) − R(j)| (here i and j are the monomer indices of two points along the primitive path). Evidently r(i,j) vanishes for i = j and is symmetric under interchange of i and j. We simulate ring polymers, so the function r(i,j) is periodic in its arguments. As we simulate in 3d periodic boundary conditions, we take care to use the nearest point among all periodic images when computing the distances |R(i) − R(j)|. The ICE primitive path has a favorable property for defining persistent contacts in that it has a noticeable persistence length (see Figure 1), which makes contacts arising from simple unentangled self-loops rare. Larger unknotted configurations such as would result from a long “hernia” (hairpin configuration of one chain segment in a melt of others) are likewise entropically unlikely to form. Because the primitive path is quite smooth on the scale of a few monomers, to speed up later calculations we slightly coarsegrain the primitive path, by partitioning the original path of N = 800 points into a sequence of 160 sets of five points each and creating a new path of 160 steps in which each point is the center of mass of five points on the original path. Thus, we represent the distance map of a given primitive path as a 160 by 160 symmetric square matrix of distance values, for which it is sufficient to display the upper triangle.
Figure 2. Folded isoconfigurational primitive paths based on different averaging times: (a) 0.1τe, (b) 0.2τe, (c) 0.4τe, and (d) 0.8τe.
for entanglement points by points where the primitive path nearly intersects itself. Note that since under Rouse scaling we have τ ∼ N2, we may say that for a choice of τa = 0.2τe a strand of (0.2)1/2Ne ≈ 0.45Ne monomers will relax during τa That is, what seems like a much shorter time than τe in fact corresponds to a modestly shorter averaged strand than Ne.
■
CLOSE CONTACTS ON THE PRIMITIVE PATH Chain-shrinking methods reduce an entangled melt to a set of piecewise linear chain paths connecting binary contact points at which two uncrossable chain segments are entangled. Thus, entanglement points are naturally identified, at the expense of reducing the full ensemble of fluctuating trajectories to a set of motionless line segments. To identify entanglement points without employing chainshrinking methods, we need a heuristic for what constitutes an
Figure 3. Distance map for a given primitive path with τa = 0.2τe plotted as a height function (a) and as a contour plot (b). E
dx.doi.org/10.1021/ma5010315 | Macromolecules XXXX, XXX, XXX−XXX
Macromolecules
Article
Figure 4. (a) Basin map corresponding to contact map of Figure 3. (b) 2d plot of basins, colored according to depth (red = deep, blue = shallow).
Figure 3a shows the distance map for a given primitive path from our simulations, which has the appearance of a rugged mountainous landscape, with many valleys separated by the mountain ridges. Evidently, the distance map has many local minima, some deeper than others. Figure 3b shows the corresponding contour plot of the distance map of Figure 3a, in which dark regions represent small distances. The contour plot clearly shows that the valleys in the distance map have a typical size in the 2d plane, i.e., a characteristic width in number of monomer points along the primitive path. Consider the pair of strands colored red and green in Figure 2, which approach close to each other and evidently form an entanglement point. The corresponding region on a distance map will be a deep valley, in which the minimum of the valley will correspond to the closest approach between the two strands. However, whereas all close entanglement points should correspond to valleys in the distance map, not all valleys in the distance map correspond to close entanglements. Two curved strands on the primitive path far from each other may also give rise to a valley on the distance map as the strands veer toward each other and then veer away again. What we expect is that deep valleys in the distance map correspond to close approaches between two strands on the primitive path. The above argument motivates us to identify all the valleys in the distance map and label them by their minimum depth. For each valley, there is one minimum point, and for each minimum point, there is a basin of attraction, such that a point moving always downhill in the distance map will ultimately reach the minimum point. We identify the set of basins of attraction for a distance map by using the conjugate gradient method to find the minimum corresponding to each starting point in the 2d grid of points in the i, j plane. The entire i, j plane is thus partitioned into a set of basins, with each basin identified by its corresponding minimum. Having constructed the basins and their corresponding minima, we can replace the distance map with a piecewise constant “basin map” in which each point in the i, j plane is replaced by the value of the minimum to which it flows. Figure 4a displays the basin map based in Figure 3. The basin map has the appearance of a landscape composed of many “basalt columns”, like Giant’s Causeway (on the coast of Northern Ireland, in County Antrim). In Figure 4b, the basin map is
viewed from above, with each basin colored according to the height of its minimum. In this view, we see how the distance map is partitioned into basins; the boundaries separating adjacent basins are the “mountain ridges” in the original distance map of Figure 3a. As previously stated, we expect that deep minima of the distance map correspond to entanglement points. We argue further that the basin of attraction for such a minimum indicates the domain of pairs of points along the primitive path that may be regarded as part of the entanglement point. Starting at one of these pairs of points, and moving along both primitive path segments at once to decrease the distance between the segments, we arrive at the corresponding minimum. Transforming the continuous distance map to the discretely valued basin map effectively discretizes the distance map and its domain, which allows us to identify entanglement points with low-lying basins. Of course, not all minima in the distance map and hence not all basins correspond to entanglement points. We must specify some criterion to determine which basins are associated with entanglement points. The simplest such criterion is to specify a maximum value r* for the depth of the minima, such that all deeper minima in the distance map are taken to be actual entanglement points. If we take r* too small, we will accept very few basins as entanglement points and miss many entanglement points that should have been included. If we take r* too large, we will accept too many basins as entanglement points. To decide what value of r* to choose, we examine the average number of basins per primitive path point M(r) with minimum depth less than r. We obtain M(r) as follows. First, we generate coarse-grained primitive paths for 100 different starting configurations from our simulation of a single N = 800 ring and construct the basin map for each. Then, for each point i on the coarse-grained primitive path, we find the number of basins with minimum depth less than r. We average this value over i and over the different primitive paths to determine M(r). (In counting the number of basins, we reject any basin for which the minimum is on the line i = j, which corresponds not to a close approach of distinct primitive path segments, but to the trivial closeness of a single primitive path segment to itself.) F
dx.doi.org/10.1021/ma5010315 | Macromolecules XXXX, XXX, XXX−XXX
Macromolecules
Article
Figure 5 displays our results for M(r), which monotonically increases from zero at r = 0 up to a plateau value of about M(r)
Figure 6. Number of local minima Nmin versus cutoff distance rc (black symbols) and linear fit (dashed line).
1.4σ, which gives both a reasonable average number of entanglement points over the entire primitive path and the pleasing result that a typical point on the primitive path participates in about two entanglement points. By properly choosing the cutoff distance rc, we can ensure that our method gives on average the correct number of entanglement points (N/Ne) over the entire primitive path. Of course, we would also like to see where on the primitive path these entanglement points are located. Figure 7 shows in (a) the set of basins corresponding to entanglement points and in (b) the complete primitive path together with periodic images of strands that form entanglement points.
Figure 5. Average number of basins M(r) versus distance r.
= 10 at about r = 5σ. (The slight upturn of the curve beyond r = 6.5 is probably due to a combination of poor statistics and the finite size of the simulation box.) To understand the shape of M(r), we argue as follows. The primitive path is smooth on the scale of the tube diameter a and can be thought of as a sequence of N/Ne entanglement strands, which are curved segments of length Ne or so and size of about a. If we take the threshold value r to be of order half the system size, then a point on a given entanglement strand will participate in a local distance minimum with every other entanglement strand. (Most of these minima, of course, will not correspond to entanglement points.) Thus, we expect a plateau value of N/Ne − 1 in M(r) when the threshold value r increases to of order L/2. Our simulation box has L = 10.46σ, hence L/2 = 5.5 or so, and from previous work we know that Ne = 65 or so; so N/Ne should be in the range 11−13.21 By the above argument, the results of Figure 5 correspond to a value of N/Ne = 11, consistent with a value of Ne = 70 or so. As our goal is to identify close entanglement points, we expect a reasonable value of cutoff distance rc should be much smaller than rc ≈ 5σ at the plateau in M(r). For example, we might guess that on average each point on the primitive path should participate in about two entanglements, so that M(r) ≈ 2. (Intuitively, any given point on the primitive path is constrained by a neighboring strand “in front” and another neighboring strand “in back”.) This would imply rc should be chosen in the range 1.4 < rc/σ < 1.6 or so. Another way to choose a value of the cutoff distance rc for identifying entanglement points is to consider the average number of entanglement points on the entire primitive path. In our construction, this corresponds to the average number of basins Nmin(rc) with minimum depth below rc. We expect the value of Nmin to be equal to N/Ne, the number of entanglement strands along the ring, which from previous work we know should be in the range 11−13. Figure 6 shows our results for Nmin(rc) for five different values of cutoff threshold (rc = 0.6, 0.9, 1.2, 1.5, and 1.6 σ). Our results for Nmin(rc) are well fit by a straight line (dashed). Evidently, a value of rc/σ in the range 1.4−1.6 gives a value of Nmin in the range 11−12, very consistent with our prior estimates of Ne. In what follows, therefore, we have chosen rc =
Figure 7. Example of below-threshold basin map and corresponding entanglements along the primitive path. Contacts colored according to depth (red = deep, blue = shallow).
The example primitive path in Figure 7 is the same as that used in Figure 1, with an averaging time of τa = 0.2τe and a cutoff distance of rc = 1.4σ. The basins and the corresponding entanglement points on the primitive path are labeled with matching numbers. The reference primitive path (shown in the unfolded representation) is colored red; the entangling strands are colored differently for easy visualization. Each entangling strand is a periodic image of a portion of the original primitive path. Note that the entanglement points are reasonably well spaced out along the reference primitive path, consistent with an intuitive picture of the primitive path as a sequence of entanglement points every Ne or so along the chain. Also, note that the basins corresponding to entanglement points have a characteristic size in the i−j plane, the linear dimension of which corresponds to about 10% of the full range of i or j. In other words, the number of monomers actively engaged in one entanglement point is about 10% of the total. This is G
dx.doi.org/10.1021/ma5010315 | Macromolecules XXXX, XXX, XXX−XXX
Macromolecules
Article
Figure 8. Top: time evolution of basin map. Basins are colored according to depth (red = deep, blue = shallow). Bottom: basins with minimum depth less than rc = 1.4σ.
configuration and perform a conventional MD simulation (with chain crossing turned off). We sample the resulting trajectory at various time intervals Δt from the starting point to obtain initial configurations for use in an ICE determination of the primitive path. In this way, we obtain a sequence in time of primitive paths for the same ring melt. For each of these primitive paths, we find the distance map, deep basins, and entanglement point locations. As a first look at the persistence of entanglement points, we examine the basin map as a function of time for a given starting configuration of our ring melt. Figure 8 (top) shows the basin maps of the initial configuration, followed by the basin map for the configuration after time delays of 2τe, 4τe, and 8τe. The bottom row of images shows the deep basins with minimum depths less than the cutoff value of rc = 1.4σ, which correspond to entanglement points. Although the set of deep basins is not completely invariant in time, we certainly notice a persistent tendency for deep basins to be located in the same places at different time delays. To further visualize the persistence of deep basins corresponding to entanglement points, we create a threedimensional plot formed from a stack of deep basin maps at different time delays. Each deep basin map forms a time slice of the set of “world lines” of entanglement points (see Figure 9). In the figure, the vertical direction is the time axis, which covers a range of time delays up to 40τe. With Ne = 65 and hence Z = 12.3, the Rouse time for the ring is τR = (Z/2)2τe ≈ 38τe (for rings, the result is (Z/2)2τe, not Z2τe as for linear chains). Hence, Figure 9 covers about one Rouse time. Strongly persistent deep basins are evident as unbroken world lines. Also visible are some deep basins that “chatter” on and off, sometimes present and sometimes absent, but repeatedly appearing at nearly the same place. These basins presumably correspond to less strongly interacting primitive path strands, which may be pulled against each other only slightly by tension along the primitive path. For such points, thermal fluctuations of the tension may be sufficient to engage and disengage the close contact of the strands, so that the corresponding deep basin minimum sometimes meets the criterion rc < 1.4σ and sometimes fails to be counted. Note that the primitive paths in ref 20 were time-averaged over much
reasonable, since we expect an entanglement point to involve about Ne monomers, and our ring is about 12Ne long (estimated from our previous value Ne = 65, so that Z = 800/65 = 12.3).21 In Figure 7, we have illustrated entanglement points as pairwise events in which only two primitive path strands participate, by selecting the appropriate strand periodic image that closely approaches the reference primitive path. However, it occurs with some small probability that three or more strands simultaneously approach each other closely. Suppose that strands near primitive path points i, j, and k are simultaneously close (with i < j < k, without loss of generality). In this case, the map of deep basins will exhibit basins at (i, k), (i, j), and (j, k), which appear roughly as a right triangle in the 2d plane (with i, k at the right angle). In Figure 7, basins 3, 4, and 10 appear in roughly this arrangement; however, for simplicity, we have omitted the third periodic image strand in each of the corresponding locations. The question of whether entanglement interactions between primitive path strands are pairwise arises as well in other recent work in our group in which we investigate primitive path motions following constraint release events.23
■
PERSISTENCE OF ENTANGLEMENT POINTS We have shown that we can state a criterion to select close approaches between two strands along the primitive path, which we identify as entanglement points. We can set the cutoff radius rc so that we find on average N/Ne such points, which do seem to correspond visually to places along the primitive path where two strands curve about each other as if under tension and restrained by uncrossability from relaxing. In addition, we have certain expectations regarding the dynamics of the tube and hence the primitive path and entanglement points. We expect that for a self-entangled ring melt, with chain crossing moves turned off, the tube should be permanent. Correspondingly, the overall shape of the primitive path should be fixed, and the location of entanglement points should persist in time. In this section, we investigate by several approaches whether the entanglement points we heuristically identify are long-lived. To do this, we start with one well-equilibrated ring H
dx.doi.org/10.1021/ma5010315 | Macromolecules XXXX, XXX, XXX−XXX
Macromolecules
Article
Figure 10. Family of deep basin autocorrelation functions g(r,t) for delay time values t/τe = 0, 0.2, 0.4, 0.8, 1.6, 3.2, and 6.4. Inset: g(0,t) versus delay time. Dashed curve is a guide to the eye.
The inset to Figure 10 plots the self-correlation g(0,t) semilog versus time, which measures the degree to which deep basins identified with entanglement points persist in time. It is apparent from the inset that this self-correlation does not decay exponentially; instead, the results suggest a persistent selfcorrelation of contact points in time, consistent with the idea of a permanent tube for our uncrossable, permanently selfentangled ring. The simple correlation function we have defined is evidently an imperfect measure of contact persistence, in that it makes no attempt to track the identity of a given contact as the participating monomers move along their respective primitive paths by Rouse motion or reptation. To do better would require recognizing the “same” persistent contact at different timesteps, which is made difficult by the phenomenon of “chattering” on and off of weak contacts.
Figure 9. Deep basin areas versus time for a given starting configuration of our ring melt, which reveal the “world lines” of entanglement points in the 2d plane of primitive path monomer indices.
longer than the ring Rouse time. The “chattering” phenomenon would not be expected to be present with such long timeaveraging, but rather reflects fluctuations in the tube on shorter time scales. Finally, we can measure the persistence of entanglement points by computing the time-dependent autocorrelation function g(r,t) for the deep basin map, as follows. We first define the characteristic function χ(r,t) for the deep basin map, such that χ(r,t) equals unity if the point r in the 2d plane of monomer indices lies within a deep basin and is zero otherwise. Then, we define the autocorrelation function as g (r, t ) = ⟨χ (0, 0)χ (r, t )⟩ − ⟨χ ⟩2
■
TIME EVOLUTION OF PRIMITIVE PATHS For a self-entangled ring polymer melt, in which the chain does not cross itself, the topological state of the system is permanent. Hence, the overall shape of the tube, which is a consequence of uncrossability, and the primitive path, which is the centerline of the tube, should be unchanging in time. However, our operative ICE definition of the primitive path is quite local in time, averaged only over a time τa of order τe. In the time regime τe < t < τR, contour length fluctuations within the tube lead to tension fluctuations along the primitive path. We expect that these fluctuations move the primitive path to some extent, as the entangled chain strands respond to the varying tension. Figure 11 illustrates how the primitive path evolves in time. Each row of figures corresponds to a different choice of averaging time τa; within each row, the primitive path of the same initial configuration at t = 0 (in red) is compared to the primitive path for the ring configuration after a time Δt has passed. (The averaging time τa increases from the top row to the bottom; the delay time Δt increases from left to right.) Evidently, longer averaging time leads to a smoother primitive path which changes less with time. In the figure, the essential shape of the primitive path varies only modestly with delay time t over the range of values considered, which extends up to about τR/3. Qualitatively, the magnitude of fluctuations in our ICE primitive path with delay time t appear somewhat smaller than
(5)
In eq 5, the angle brackets ⟨ ⟩ denote averaging over the choice of 2d origin (i.e., the starting location) and over the temporal origin (i.e., the starting time). Because the 2d domain of primitive path indices is periodic in both directions, in practice we can use discrete Fourier transform methods to evaluate the spatial correlation. (We first apply a discrete Fourier transform to χ(r,0) and χ(r,t) to obtain χ(k,0) and χ(k,t); the desired correlation is then the inverse discrete Fourier transform of the product χ*(k,0)χ(k,t).) Although g(r,t) is actually a function of the two-dimensional displacement r between two points in the 2d plane of monomer indices, we compute an isotropic average g(|r|,t) ≡ g(r,t) by averaging over all displacements r with the same length. Figure 10 displays a family of g(r,t) curves for time delays t equal to 0, 0.2, 0.4, 0.8, 1.6, 3.2, and 6.4 times τe. It is apparent from the figure that the locations of the deep basins are only weakly spatially correlated (there is only a modest wiggle in the correlation function indicating near-neighbor correlations); however, the correlation of a basin with itself (peak at r = 0) does persist in time. I
dx.doi.org/10.1021/ma5010315 | Macromolecules XXXX, XXX, XXX−XXX
Macromolecules
Article
Figure 11. Time evolution of primitive path from its initial shape (red curve) for different delay times t and different averaging times τa.
dynamics of the p = 0 Rouse mode. We remove the contribution of ring reptation by defining the distance between two primitive paths as the minimum of the average square displacement
one might expect from discussion in ref 20. One possible origin for the difference is that we are simulating a single self-entangled ring of N = 800 in 3d periodic boundary conditions, rather than a system of hundreds of rings of comparable size. As a result, our system is much smaller, and the cutoff wavelength for “gel modes” is shorter. Figure 11 shows that the path does vary in time, apparently around a well-defined permanent shape, presumably in response to tension fluctuations. To quantify the size of displacements of the primitive path with time and to extend the range of time delays over which we measure such displacements, we gather statistics on the displacements of the primitive path, described by the location of successive points Ri(t). In quantifying the distance between two primitive paths separated in time, we must account for the fact that the ring polymer can move within the tube and hence move along the primitive path. Such motions change the monomer labels along the path, but not the shape of the path itself. We may distinguish two sorts of motion of the ring along the primitive path: curvilinear Rouse motion, in which the ring “bunches up” and “stretches out” along the path, and reptation, in which the ring moves “rigidly” along the path, i.e., without bunching up or stretching out. Curvilinear Rouse motion has the effect of changing the spacing of the primitive path points along the path, while reptation has the effect of uniformly shifting the monomer labels along the path. Both motions leave the path itself undisturbed but change the function Ri(t) we use to describe the path. It is awkward to remove the effect of the curvilinear Rouse motion on the apparent displacement of the primitive path. (This is accomplished by the complete arc length remapping in the tube axis algorithm of ref 20, which we have not employed here.) However, it is straightforward to remove the effect of reptation, which for rings corresponds simply to the diffusive
rs 2(t , δN ) =
1 N
N
∑ (R i+ δN(t0 + t ) − R i(t0))2 i=1
(6)
with respect to the shift δN of the (periodic) monomer indices. This simple “uniform shift” remapping may be regarded as a “poor man’s” version of complete arclength remapping between the two primitive paths. Figure 12 illustrates this calculation, by plotting the distance between two primitive paths as a function of δN, for example pairs of primitive paths separated in time by t with t/τe = 0.2, 1,
Figure 12. Difference between two primitive paths from eq 6 versus monomer index shift δN for different time intervals Δt (□: 0.2τe; ○: τe; △: 4τe; ▽: 20τe). J
dx.doi.org/10.1021/ma5010315 | Macromolecules XXXX, XXX, XXX−XXX
Macromolecules
Article
4, and 20. The figure shows how rs2(t,δN) displays a clear minimum, at a value of δN which is quite small for short time intervals but tends to be larger for longer time delays, as expected. By averaging the quantity minδN rs2(t,δN) over many instances of primitive paths separated in time by t, we obtain the results of Figure 13. The range of delay times considered
urations by adding to our molecular dynamics simulations various Monte Carlo moves that allow the chain to cross itself. Our system is ideal for studying the tube and its properties because entanglements are permanent for rings, once chain crossing moves are turned off. (We study one ring rather than a melt of many rings, to minimize any corrections that arise from the finite length of the rings.) In identifying points along the primitive path where strands closely approach each other, we take advantage of the ability to choose the isoconfigurational averaging time τa. It turns out that by choosing τa = 0.2τe, our averaging procedure leads to tight entanglement points appearing as near intersections between primitive path strands. We then construct the map of squared distances δr2(i,j) between primitive path points i and j, find the basins of attraction for each local minimum, and select among these basins those with minimum distances closer than a threshold value rc to identify points of close approach. We select the distance threshold rc = 1.4σ, so that on average we find the expected number N/Ne of entanglement points, equal to the number of entanglement strands in our system. (From previous work, we have a good value for Ne, around Ne = 65.) For this choice of rc, we find also that a typical point i0 on the primitive path is within the basin of attraction of about two entanglement points. To be within the basin of attraction means: if we move continuously along the primitive path starting from i0 to some nearby i, we will be moving closer in space to some corresponding closest approach point j on another portion of the primitive path. By adjusting the threshold distance rc, we can obtain the expected number of entanglement points. A test of whether these points are meaningful is to see whether they persist in time as the conformation of our self-entangled ring fluctuates in time. To test this, we generated a sequence of primitive paths of a self-entangled ring melt as a function of time delay from a single starting configuration. For each of these primitive paths, we determined a set of deep basins and corresponding locations of entanglement points. By stacking the two-dimensional plots of the deep basin areas, we generate a 3d representation of the “world lines” of the entanglement points. This visualization makes clear that the deep basins corresponding to entanglement points persist in time to some extent for our self-entangled ring, although entanglement points corresponding to some more weakly interacting strands do tend to “chatter” on and off, as the tension along the strands fluctuates. Indeed, any definition of persistent contacts via a threshold in distance or force will be subject to such “chattering”, for contacts fluctuating about the threshold. This will be true regardless of how long we time-average or otherwise smooth the primitive path, whereas the primitive path itself defined by ICE averaging or similar methods will have a continuous dynamical trajectory. In this sense persistent contacts would seem to give a less robust description of entangled chains than the time-averaged primitive path, however conceptually convenient or appealing a set of discrete pairwise entanglement points may be. We also compute the time-dependent correlation function g(r,t) of the deep basin characteristic function χ(r,t) (χ(r,t) = 1 if a pair of path points r = {i, j} is in a deep basin at time t and 0 otherwise). The autocorrelation function quantifies the persistence of the deep basin mapwith the significant limitation that as defined, it does not account for collective motion of the chain inside its tube, whether by reptation or Rouse motion, that relabels the monomers participating in a
Figure 13. Averaged minimum square distance ⟨minδN rs2(t,δN)⟩ between primitive paths separated in time by t.
extends to 106τLJ, comparable to the duration of simulations reported in ref 20. The average square distance between paths grows until a plateau is reached, at a time around 105τLJ. We can identify the time scale at which the primitive path displacement stops growing, by comparing to our results for the mean-square displacement of ring monomers. All monomers on a ring are statistically equivalent, so it suffices to obtain ⟨ΔR2(t)⟩ averaged over all monomers. With our value of Ne ≈ 65 and hence Z ≈ 12, we can then compute the Rouse time for the ring of τR = (Z/2)2τe (for rings, the result is (Z/2)2τe, not Z2τe as for linear chains). With our value of τe of about 2000τLJ,24 this gives τR ≈ 8 × 104τLJ, consistent with the time scale on which the primitive path appears to stop changing. Since our primitive path has not been time-averaged beyond τe, and we have not implemented the full remapping procedure of the tube axis algorithm, it is reasonable that there should be fluctuations in the primitive path position up to delay times of order τR, either because of physical fluctuations of the tube or because of fluctuations in monomer positions along the tube.
■
CONCLUSIONS There are two complementary (sometimes competing) views of entanglement effects in polymer melts and solutions: chains confined to a continuous tube or constrained by discrete entanglement points. In this paper, we have presented a heuristic approach for finding entanglement points in simulations of entangled polymer melts by looking for close approaches between two strands along the primitive path. In our work, we identify the primitive path by means of an isoconfigurational ensemble (ICE) averaging of many short simulation trajectories starting from the same chain configuration, but with different random initial velocities. In this way, the ICE primitive path averages over all the ways the chain might wiggle within its confining tube. We study melts of one long self-entangled ring polymer in periodic boundary conditions, as a proxy for entangled melts of long linear chains. We equilibrate the entanglement configK
dx.doi.org/10.1021/ma5010315 | Macromolecules XXXX, XXX, XXX−XXX
Macromolecules
Article
(24) Qin, J.; Milner, S. T., to be published. (25) Likhtman, A. E.; Ponmurugan, M. Macromolecules 2014, 47, 1470−1481. (26) Qin, J.; Morse, D. C. J. Chem. Phys. 2009, 130, 224902.
given contact. In future work, we can improve on our definition by maximizing the correlation function at each time with respect to cyclic relabeling of the monomers, which will compensate for reptative displacements within the tube. We are likewise interested in the persistence in time of the primitive path itself. Since the topological state of our selfentangled ring melt is permanent once chain crossing moves are turned off, we expect the overall shape of the tube to be likewise permanent. However, our ICE definition of the primitive path is local in time, in the sense that we average over rather short simulations of duration τa of order τe to identify the path. In the time regime τe < t < τR, contour-length fluctuations within the tube are expected to shuffle monomers around in the tube and give rise to tension fluctuations along the primitive path. Ultimately, these contour-length fluctuations give rise to fluctuations in the primitive path location, which appear to scale with Rouse-like dynamical exponent ΔR2(t) ∼ t1/2 in the regime τe < t < τR, with mean-square displacements of order the tube diameter at the end of the regime.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected] (S.T.M.). Notes
The authors declare no competing financial interest.
■ ■
ACKNOWLEDGMENTS The authors acknowledge support from NSF DMR-0907376 for this work. REFERENCES
(1) de Gennes, P. G. J. Chem. Phys. 1971, 55, 572. (2) Doi, M.; Edwards, S. The Theory of Polymer Dynamics; Oxford University Press: New York, 1988. (3) de Gennes, P. G. J. Phys. (Paris) 1975, 36, 1199. (4) Rubinstein, M.; Colby, R. H. J. Chem. Phys. 1988, 89, 5291−5306. (5) Milner, S. T.; McLeish, T. C. B. Phys. Rev. Lett. 1998, 81, 725− 728. (6) Cloizeaux, J. D. Macromolecules 1990, 23, 4678−4687. (7) Likhtman, A. E.; McLeish, T. C. B. Macromolecules 2002, 35, 6332−6343. (8) McLeish, T. C. B. Adv. Phys. 2002, 51, 1379−1527. (9) Milner, S. T.; McLeish, T. C. B.; Likhtman, A. E. J. Rheol. 2001, 45, 539−563. (10) Graham, R. S.; Likhtman, A. E.; McLeish, T. C. B.; Milner, S. T. J. Rheol. 2003, 47, 1171. (11) Likhtman, A. E.; Graham, R. S. J. Non-Newtonian Fluid Mech. 2003, 114, 1−12. (12) Lin, Y. H. Macromolecules 1987, 20, 3080−3083. (13) Milner, S. Macromolecules 2005, 38, 4929−4939. (14) Everaers, R.; Sukumaran, S. K.; Grest, G. S.; Svaneborg, C.; Sivasubramanian, A.; Kremer, K. Science 2004, 303, 823−826. (15) Zhou, Q.; Larson, R. G. Macromolecules 2005, 38, 5761−5765. (16) Kröger, M. Comput. Phys. Commun. 2005, 168, 209−232. (17) Tzoumanekas, C.; Theodorou, D. N. Macromolecules 2006, 39, 4592−4604. (18) Read, D. J.; Jagannathan, K.; Likhtman, A. E. Macromolecules 2008, 41, 6843−6853. (19) Bisbee, W.; Qin, J.; Milner, S. T. Macromolecules 2011, 44, 8972−8980. (20) Likhtman, A. E. Soft Matter 2014, 10, 1895−1904. (21) Qin, J.; Milner, S. T. Soft Matter 2011, 7, 10676. (22) Qin, J.; Milner, S. T. Macromolecules 2014, 47, 6077−6085. (23) Cao, J.; Qin, J.; Milner, S. T. Macromolecules 2014, 47, 2479− 2486. L
dx.doi.org/10.1021/ma5010315 | Macromolecules XXXX, XXX, XXX−XXX