Tube Diameter of Oriented and Stretched Polymer Melts

Article pubs.acs.org/Macromolecules

Tube Diameter of Oriented and Stretched Polymer Melts Jian Qin* and Scott T. Milner Department of Chemical Engineering, The Pennsylvania State University, University Park, Pennsylvania 16802, United States ABSTRACT: The Lin−Noolandi entanglement ansatz correlates the tube diameter with polymer density and chain stiffness, and has been successfully extended to concentrated solutions as well. We develop a new version of the entanglement ansatz, applicable to polymer melts under uniform tension, and find that the tube diameter equals its unperturbed value a when the pulling force F is less than the thermal tension kT/a, and scales as F−1/3 when F exceeds kT/a. We verify this behavior by using isoconfigurational ensemble averaging to estimate the tube diameter in simulated polymer melts that are oriented by pulling on the ends of topologically equilibrated chains. When the tension is large enough to nearly fully extend the chain, our analysis predicts that the tube diameter scales as F−2/3, but our simulation results do not reach this regime.

■

INTRODUCTION The tube concept lies at the heart of modern theory of polymer rheology.1 It originates from the fact that two chains cannot cut through each other in the course of motion. This molecular uncrossability effectively confines the polymer chain into a tube-like region, and severely restricts its motion along the directions orthogonal to the tube contour. The strength of this effect is measured by the tube diameter a, or equivalently the entanglement length Ne, which is related to a by random walk statistics (a = Ne1/2b, where b is the statistical segment length). The tube diameter is one of the two input parameters required by rheological theories.2 In practice, values of a are obtained by comparing predictions of rheological theory and experimental measurements, particularly the plateau modulus for long entangled linear melts.3 The tube diameter is an emergent length scale, dependent on microscopic material parameters such as monomer bulkiness and chain stiffness, but not on chain length or molecular topology (such as long chain branching). Currently, no mechanistic theory exists that predicts the tube diameter from microscopic material parameters. However, an empirical relation based on the Lin−Noolandi entanglement ansatz4,5 successfully correlates the tube diameter with chain dimension and displaced volume, for a wide range of entangled polymer melts.3 With a different ansatz, Colby and Rubinstein developed a scaling relation that described the concentration dependence of the tube diameter for concentrated solutions.6 These two approaches were reconciled in ref 7, which developed a single entanglement ansatz that predicted how the tube diameter depends on both chain concentration and chain architecture (monomer bulkiness and chain stiffness). Practically, this relation is sufficient for calculating the flow and deformation behaviors of polymers in the linear regime. However, for flows and deformations with large amplitudes, the polymer chains are oriented at the molecular scale. This perturbation at the molecular scale will potentially change the tube diameter from its value in the unoriented isotropic state. Any tube-based theory of nonlinear rheology must make some © 2013 American Chemical Society

assumption about how the tube diameter varies with chain orientation and stretching. The GLaMM model,8 a highly successful theory of polymer rheology derived from the tube model, makes an explicit assumption about how the tube diameter changes in fast flows that both the entanglement length and tube diameter decrease with chain stretching. However, this assumption has not been directly validated. Neither has this assumption been shown to be consistent with the entanglement ansatz of ref 7, which led to a unified account of the tube diameter in isotropic melts and solutions. This is an important gap in our phenomenological understanding of the tube diameter. In this work, we study how the tube diameter varies with chain stretching. The work consists of two parts. First, we develop three competing versions of entanglement ansatz for polymer melts consisting of chains under tension. Two of these are straightforward applications of the Lin−Noolandi and Colby−Rubinstein assumptions, and the third is a novel ansatz. The three approaches all agree for isotropic melts, but differ in their predictions for diluted melts and melts of chains under tension. The results are presented in the second and third sections. To investigate which (if any) of these proposals are correct, we conduct Monte Carlo simulations of oriented polymer melts, use isoconfigurational ensemble averaging to estimate the tube diameter,9 and compare the results with our predictions. The details of our model and simulation algorithms are summarized in the fourth section, and the simulation results are presented in the fifth section. The results section is divided into five subsections: the first discusses the force−extension relation for chains in melts; the second presents the main simulation result on the force dependence of tube diameter; the third presents an extension of the theory of the third section to account for the simulation Received: October 5, 2012 Revised: December 3, 2012 Published: February 13, 2013 1659

dx.doi.org/10.1021/ma302095k | Macromolecules 2013, 46, 1659−1672

Macromolecules

Article

melts,7 and to stretched melts, as we will see below. In ref7., it was shown that ansatz I implies the entanglement molecular weight in diluted melts would scale with polymer volume fraction ϕ as ϕ−2, whereas using ansatz II results in a scaling of ϕ−4/3. (See Appendix for a brief reprise of these scaling results.) Ref 7 chose ansatz II since the ϕ−4/3 scaling agrees with earlier experimental results in semidilute Θ solutions.6,10 As a consequence, the plateau modulus G0(ϕ) for semidilute Θ solutions (such as result from mixtures of long entangled and short unentangled chains) scales as ϕ7/3, which seemed to be consistent with available rheological data. In contrast, ansatz I implies that G0(ϕ) scales as ϕ3, which is much too strong a dilution effect when compared to available data. However, subsequent experiments have cast doubt on the validity of the ansatz II scaling; a recent thorough discussion by Larson et al.11 suggests that Ne(ϕ) appears to scale linearly with 1/ϕ (corresponding to G0(ϕ) scaling as ϕ2). This scaling would result if we regard entanglements as localized binary contacts along the length of a given chain, and the entanglement molecular weight as the arclength between successive entanglements. In a diluted melt, only a fraction ϕ of the contacts would be present (or “operative”, in a mixture of long entangled chains and short quickly relaxing chains), hence Ne(ϕ) would scale as 1/ϕ. More recently, we have carried out simulations to investigate the dilation effects on the tube diameter, obtained both by direct observation, and by topological methods;12 our recent results also support the ϕ−1 scaling for the entanglement molecular weight. To reconcile this scaling with success of the original Lin−Noolandi idea as applied to isotropic melts, we introduce ansatz III, in which we assert that an entanglement strand consists of a fixed number of close contacts with other chains, and each close contact is the size of a packing blob (volume p3, with Np monomers). We express this as

data; the fourth compares our results to previously assumed tube diameter scalings in the rheology literature; the fifth analyzes the behavior of the tube diameter in the strong stretching limit. We summarize and discuss our findings in the last section.

■

ENTANGLEMENT ANSÄ TZE We first summarize the original Lin−Noolandi ansatz. In a dense polymer melt, each strand of monomers cohabits a volume of space along with many strands of similar length. The longer a given strand, the more neighboring strands cohabit the same pervaded volume. The ansatz asserts that when the number of neighbors cohabiting the same pervaded volume reaches a constant C, an entanglement results. The corresponding strand length is identified as Ne, with the tube diameter a identified as the mean-square end-to-end distance of an entanglement strand, a2 = Neb2.7 This statement translates to the following C=

a3 Ne Ω

(1)

in which a 3 estimates the volume pervaded by one entanglement strand and Ω denotes the monomer volume. Then the tube diameter can be written as

a=C

Ω = Cp b2

(2)

in which we have introduced the packing length p defined as p = NΩ/R2 = Ω/b2, where N is the chain length and R the chain end to end separation vector. The value of p increases for bulky chains (large Ω) and decreases for stiff chains (large b2). The tube diameter is proportional to the packing length. Correspondingly, the entanglement length can be written as ⎛ Ω ⎞2 Ne = C 2⎜ 3 ⎟ = C 2Np ⎝b ⎠

ϕNe(ϕ) = C2 Np

(3)

in which we have introduced the number of monomers per “packing blob” Np, defined as NpΩ = p3. The packing length is the length scale below which a random walk overfills space,7 i.e., it is the length scale of a strand of Np monomers such that (Np1/2b)3 = NpΩ. Below this length scale, the space around a given monomer is primarily filled with monomers that are on the same chain and are close along the arc length. Thus, the close contact of the given monomer with other monomers is excluded below the packing length. This suggests that the volume of a close contact between polymer strands may be estimated as p3. Several different versions of an entanglement ansatz can be formulated, all of which are equivalent for melts of undeformed chains. The first (Lin−Noolandi) version, where the pervaded volume of an entanglement strand contains a fixed number of neighboring strands, we denote as ansatz I. Alternately, we may assert (following Colby and Rubinstein6 and ref 7) that the pervaded volume of an entanglement strand contains a fixed number of close contacts, namely a 3 / p3 = C 3

(5)

which reduces to eq 3 for ϕ = 1. In the next section, we apply all three versions of the entanglement ansatz to chains in a melt under tension, and determine how the Ne scales with chain tension for each version.

■

ENTANGLEMENT ANSÄ TZE FOR CHAINS UNDER TENSION To estimate the entanglement length for stretched chains, we first recall how chain conformations are perturbed by tension. A chain under tension F can be regarded as a sequence of tension blobs13 of dimension ξ. The chain conformation is approximately a linear string of blobs, aligned with the direction of force. More precisely, the sequence of blobs adopts randomwalk configurations in the plane normal to the force. Each blob is an approximately unperturbed random walk of g monomers, satisfying ξ2 = gb2. The relation between the applied force F and the blob size ξ can be written as

(4)

using eq 2. We denote this version as ansatz II. Ansatz I and ansatz II are equivalent for undeformed polymer melts, but yield different results for the tube diameter and entanglement molecular weight when generalized to diluted

F=

3kT ξ

(6)

The end-to-end distance becomes 1660


Macromolecules L = (N / g )ξ =

Article

Nb2 F = NF /K 3kT

Now we determine the tube diameter, which sets the scale of transverse fluctuations of a chain in its tube. In an isotropic melt, the tube diameter is taken to be the mean-square end-toend distance of an entanglement strand, and the pervaded volume of an entanglement strand scales as a3 . For entanglement strands under tension, the tube diameter sets the scale of transverse fluctuations of an entanglement strand, which we identify with the cylinder radius D = n1/2ξ; this leads to a a(f ) = (aξ)1/2 = 1/2 f (11)

(7)

2

in which K = 3kT/b is the correct spring constant per monomer, validating the relation between F and ξ. The larger a force is applied, the smaller ξ becomes. The chain is linearly extended at length scales above ξ, and is largely unperturbed at smaller length scales. For a force equal to the thermal tension FT = 3kT/a,1 we have ξ = a, and the chain begins to be perturbed at the length scale of the tube itself. Evidently we may write

a F = =f ξ FT

so that the tube diameter for chains under tension decreases, scaling as f−1/2. (The length of the cylinder, nξ = a, remains constant.) Ansatz II. Now assume ansatz II is correct, so that an entanglement is defined by a fixed number C3 of close contacts in the pervaded volume of an entanglement strand. This leads to

(8)

in which we have defined f as the applied force in units of the thermal tension. If the applied force is smaller than FT, then ξ is larger than a, and chain conformations at the scale of the tube are little affected by the applied force. Hence for weak applied forces we do not expect the statistics of entanglement to change, and the tube diameter will be little affected. If the applied force F exceeds the thermal tension FT (hence f > 1), the tension blob size ξ becomes smaller than the equilibrium tube diameter. Then chain conformations at the scale of the tube are perturbed, and we expect the various forms of the Lin−Noolandi ansatz are affected. For completeness, we consider in turn ansatz I, II, and III, even though we have argued that only ansatz III appears to be consistent with rheological and simulation results for melts and Θ solutions. Let n be the number of tension blobs needed to create an entanglement. The tension blobs execute a random walk in the transverse direction as they march along in the force direction. The pervaded volume of n blobs thus has the shape of a cylinder, with length L = nξ and diameter D = (nξ2)1/2 as shown in Figure 1.

C3 =

⎛ ξ ⎞1/2 Ne Ne(f ) = ng = Ne⎜ ⎟ = 1/2 ⎝a⎠ f

As we shall show in the results section, this dependence is too weak to capture the trend of simulation results; it is also likely too weak to give good agreement with rheological data if used in a GLaMM-type model. Ansatz III. Thus, far, we have ansatz I, which gives too strong a dependence of Ne(ϕ) for diluted melts but reasonable scaling of Ne(F) for melts under tension, and ansatz II, which gives reasonable results for diluted melts but too weak a dependence of Ne(F) for melts under tension. Of course, we want a single approach to describing entanglement scaling that works for melts, solutions, and melts under tension. As expressed above, ansatz III assumes that an entanglement strand is defined by a fixed number of close contacts with neighboring strands. The close contacts are assumed to have a size of order the packing volume p3. We are not applying a force nearly large enough to perturb chain conformations on the scale of p. We can regard our chains under tension as biased random walks of unperturbed packing blobs. Hence a naive interpretation of ansatz III would imply that an entanglement strand under tension would consist of the same number of packing blobs, and thus that Ne would be unaffected by the tension. Instead, we reconsider ansatz III, from the point of view of topological complexity. In a melt under tension, the chain trajectories are biased random walks. The chain trajectories consist of an average linear displacement with increasing s along the z axis, with transverse and longitudinal fluctuations superimposed. If we imagine traveling along the arclength of

(9)

Using the relations a = Cp and f = a/ξ, we find that the number of tension blobs in an entanglement strand n equals f. Then the entanglement strand length is Ne(f ) = ng = Ne

N ξ = e a f

(13)

The corresponding scaling prediction for the tube diameter is again D = n1/2ξ, which this time gives a a(f ) = 1/4 f (14)

Ansatz I. Applying ansatz I, we assert that the number of strands of n tension blobs in the cylindrically shaped pervaded volume needed to create an entanglement is the same as for an unstretched melt. This leads to nξ·nξ 2 nξ = ng Ω p

(12)

Recalling that for undeformed melts we have C3 = a3/p3, and replacing ξ using f = a/ξ, we find n = f 3/2. Then the entanglement strand length is

Figure 1. Pervaded volume of a string of blobs under tension.

C=

n 2ξ 3 p3

(10)

Thus, ansatz I predicts that the entanglement length scales as f−1. Qualitatively, stretching chains out brings more of them into close proximity, rather like the effect of stiffening the chains. 1661


Macromolecules

■

a chain, the state of entanglement of a given chain with its neighbors is not altered by the average uniform motion along z of the given chain together with its neighbors. This average motion simply serves to stretch out the trajectories. We envision removing this uniform motion, by shrinking the end-to-end displacement of the chain by some factor λ; we choose λ in such a way that an entanglement strand (the length of which we do not yet know) will be roughly isotropic after the deformation. In this way, we hope to recover an ensemble of roughly isotropic random-walk configurations, for which the statistics of entanglement will be similar to an isotropic melt. We do this because we do not otherwise know how to relate the entanglement statistics of a melt of biased random walks to the entanglement statistics of an isotropic melt. Note that in shrinking the end-to-end displacement of an entanglement strand by a factor λ, we have decreased the pervaded volume of the strand by the same factor. However, we want to preserve the topological complexity of the melt under this shrinking, neither increasing nor decreasing the number of distinct knots in which the melt may be found. We assert that on shrinking the pervaded volume, we must likewise shrink the packing volume p3 by a factor of λ. Otherwise, multiple distinct topological states in the melt under tension must map to the same knotted state in the shrunken melt, because there is not enough “spatial resolution” in the shrunken melt to always tell whether a given strand passed on one side or the other of a nearby strand. In short, we make the following scaling assertions:

Us = F(R z+ − R z−) ≡ FR z

(15)

in which we have defined Np′Ω = p′3. From the first relation above, we find n = λ2. Dividing the third relation by its undeformed limit, and using the second relation to eliminate p and p′, we obtain ⎛ N (f ) ⎞−1/2 Ne(f ) = λ−1 = n−1/2 = ⎜ e ⎟ Ne ⎝ g ⎠ ⎛ Ne(f ) ⎞−1/2 ⎛ Ne ⎞−1/2 =⎜ ⎜ ⎟ ⎟ ⎝ Ne ⎠ ⎝g⎠

(16)

Rearranging, we have −1/3 ⎛ a ⎞−2/3 Ne(f ) ⎛ Ne ⎞ =⎜ ⎟ =⎜ ⎟ = f −2/3 ⎝ξ⎠ Ne ⎝g⎠

(17)

The corresponding scaling for the tube diameter, again defined by transverse fluctuations, is a(f ) = af −1/3

(19)

where F is the magnitude of forces, R+z and R−z are z components of the end bead positions. We used various Monte Carlo moves to equilibrate the system and sample configurations: hybrid MC/MD, configuration biased single-chain and double-chain rebridging moves, tethered bead displacement, and head−tail-flip move. The first three moves have been used in our previous simulations of unstretched systems.9,16 The last two moves are introduced here for simulating stretched chains. The tethered bead displacement move randomly chooses an end bead, displaces its position randomly, and accepts the proposed position with the Metropolis rule. The head−tail-flip move reverses the sign of the force on two end beads of a given chain, and is also accepted according to the Metropolis rule. This last move changes the sign of Rz of the chain, and helps to equilibrate the system. It is analogous to the chain head−tail-swap move in diblock copolymer simulations, which allows phase space to be more effectively sampled.17 Even with these moves, equilibration becomes difficult as the force value increases. To speed up the equilibration, we also implemented a version of replica exchange. In our scheme, a sequence of systems (replicas) are simulated in parallel, each with a different stretching force (typically 12−16 replicas are used). Neighboring replicas are allowed to attempt to exchange configurations periodically, accordingly to the Metropolis rule. This move effectively allows high-tension configurations to communicate with lower tension ones. All simulations were conducted using the general purpose Monte Carlo and molecular dynamics simulation package (“Simpatico”18). In simulations of unstretched systems, the initial configurations are prepared by using unbiased random walks. Before the more sophisticated Monte Carlo moves are turned on, we employ random bead displacement moves to equilibrate the system pressure. To generate seed configurations for most of the stretched systems, we begin with the unstretched system, equilibrate it, and then increase the force values progressively. At each force value, the system is given enough time to reequilibrate. To confirm that this is sufficient for generating stretched system configurations, we also used a different way to generate the seed configuration: we use the force−extension relation for Gaussian chains to estimate the extent of bond orientation, and use the Boltzmann weight to construct appropriately biased random walks. The results obtained this way are in good agreement with those obtained from progressively stretched systems.

p3 fixed topological complexity λ

Ne = C 2 fixed number of packing blobs per strand Np ′

SIMULATION MODEL AND METHODS

The key result we obtained in the previous section is that the tube diameter a scales with f−1/3 under the assumption of ansatz III (eq 18), which seems most likely given the existing data. To test this prediction, we simulated melts of entangled polymers under tension, and estimated the tension dependence of the tube diameter. We used a coarse grained bead−spring model, in which the bonded beads interact through a harmonic potential, instead of the FENE potential as in the Kremer−Grest14 model. The nonbonded beads are interacting through the purely repulsive Weeks−Chandler−Anderson (WCA) pair potential.15 To impose the uniaxial stretching, we apply a tethering force along the z direction of equal magnitudes but opposite signs to the beads at the chain ends. This corresponds to an additional potential energy in the system Hamiltonian, of the form

nξ = (nξ 2)1/2 shrink to isotropy λ p′3 =

Article

(18)

As we shall show below, a theory based on ansatz III of the crossover from weak to strong tension agrees reasonably well with our simulation results, without adjustable parameters. 1662


Macromolecules

Article

We simulated melts of M = 10 chains of N = 400 beads. For our model Ne is about 65,9 so N = 400 chains have about 6 entanglement strands. The chains fill a tetragonal box of dimension 12 × 12 × 39.7σ3, which corresponds to a bead number density ρ = 0.7σ−3, where σ is the range of the WCA (truncated L-J) pair potential. The box dimension along x and y direction is about Rg, and that along z direction is about 3.5Rg. This bead filling density has been verified to represent the melt condition properly in simulations of binary homopolymer melts.19 We used periodic boundary conditions throughout, so it is important to track bead positions in a fixed coordinate system before applying the stretch potential eq 19 in the Metropolis test. To investigate the force−extension relation in the nonlinear regime, we also simulated three comparison systems with chain length N = 50, 100, and 200 and number of chains M = 80, 40, and 20 respectively. Each of these comparison systems contain the same number of beads as the N = 400 system. The box dimensions along the z direction for these systems were taken equal to 10, 14.14, and 20 (to scale with √N), and those along the x and y directions are identical and are chosen to maintain the constant bead number density (0.7σ−3).

The oriented melts differ from the isotropic melts only at longer length scales. Individual bonds between beads are not stretched. The bond length distributions are indistinguishable in the range of force values we investigated. In contrast, the bond orientation becomes increasingly aligned along the z axis as the force increases. The bond−bond autocorrelation functions damp to zero within less than 10 bonds, even for the most strongly oriented melts we investigated. We quantify the extent of chain orientation by measuring the mean end to end distance Rz. The dependence of Rz on the applied force (force−extension law) is presented in the next subsection. In the second subsection, we describe how the isoconfigurational ensemble averaging is used to estimate the tube diameter, and present the results on the force dependence of the tube diameter. In the third subsection, we extend the three versions of the entanglement ansatz to describe the crossover between weak (F < FT) and strong (F > FT) applied tensions; we compare these crossoverexpressions to our simulation results. In the fourth subsection, we discuss the relation between our work and various previous theories. In the fifth subsection, by applying the three forms of entanglement ansatz to the freely joined chain model, we analyze the behavior of the tube diameter in the limit of full extension. Force−Extension Law. A significant amount of previous work has examined the force−extension relation for an isolated chain, in vacuum or in solution (e.g., refs 20−22). Very little work has focused on how polymers in an entangled melt respond to pulling force (Smith and Kober et al. simulated the force−extension relation for a single PDMS chain in oligomers23), presumably because of the difficulty of equilibrating oriented dense melts. Before showing our results, we first summarize the force− extension relation of a Gaussian chain. We expect chains in the melt to be reasonably well described by Gaussian random walks because of the screening of self-avoidance.1 We expect a Gaussian description to be valid as long as we do not apply too large a force, so that the chain begins to be affected by the limit of full extension. The probability distribution for the end to end separation vector R of a Gaussian chain P(R) is proportional to exp(−3kTR2/(2Nb2)). So the configuration free energy of a Gaussian chain under a uniaxial pulling force can be written as 3kTR 2 /(2Nb2 ) − FR z . Within the linear regime, the distributions of components Rx and Ry are both Gaussian with mean zero and variance Nb2/3. The distribution of Rz is also Gaussian, with mean

■

SIMULATION RESULTS The goal of our simulation work is to examine how the tube diameter changes in oriented melts, and to test the predictions of the three versions of entanglement ansatz. We generate a set of oriented melts by applying pulling forces at chain ends, as described in the previous section. (More physically, we might prefer to study a system oriented by applying a shear flow. However, a sheared system cannot be represented with an equilibrium ensemble with a well-defined Hamiltonian, which is needed by the Metropolis tests for the chain rebridging moves that we employ in our simulations.) As the applied force increases, the chain become oriented and stretched. Two typical snapshots from our simulations are shown in Figure 2; the left image depicts a system without tension, and the right shows a strongly oriented melt. For convenience of visualization, different chains are colored differently. The green and orange spheres label the chain ends pulled along the opposite directions. The two contrasting figures show clearly how chains are aligned along the z axis.

⟨R z⟩ =

Nb2 F 3kT

(20)

2

and variance Nb /3. Thus, a Gaussian chain has a strictly linear force−extension law. For small applied force, all these predictions are confirmed by our data. The distributions of Rx and Ry are described nearly perfectly by the Gaussian with the predicted variance (independently measured statistical segment length b = 1.404σ is used19). The Rz distribution also fits a Gaussian with the correct variance and nonzero mean. Figure 3 shows the distributions of Rz for our N = 400 system at different forces (dashed curves are Gaussian predictions for F = 0 and 0.115 kT/σ respectively). The results at F = 0 fits a Gaussian with variance (1/10)(Nb2/3) exactly, where the extra factor 1/10 accounts for the fact that we have 10 chains in the system. As F increases, the distribution is shifted to the right side. For F less than 0.2 kT/σ, the width of distribution is essentially constant.

Figure 2. Simulation snapshots of the N = 400 system at tensions F = 0 and 0.3kT/σ (about 1.2 times FT). Simulation cells are also shown; the long axis is the z direction. Orange and green spheres are used to label positions of chain ends pulled along the +z and −z directions, respectively. 1663


Macromolecules

Article

Figure 5. Reduced force−extension relation from results of all chain lengths. Dashed line: Gaussian prediction, Rz = Nb2F/(3kT). Black solid line: Langevin prediction, Rz = Nl03 (FlK,0/(kT)). Blue solid line: nonlinear Langevin prediction, Rz = Nl03 (FlK(F)/(kT)).

Figure 3. From left to right: distributions of Rz values at tensions F = 0, 0.115, 0.190, 0.246, and 0.277 kT/σ. For F = 0 and 0.115 kT/σ, the dashed lines are Gaussian predictions.

But for greater F values, the distribution narrows, signaling the crossover to the nonlinear deformation regime. Figure 4 shows how ⟨Rz⟩ varies with the applied force for systems with N = 50, 100, 200, and 400. The dashed lines are Gaussian predictions given by eq 20. For all N values, the simulation data follows the Gaussian prediction up to about F = 0.2 kT/σ. For larger forces, Rz values from simulation are less than predicted by eq 20. Among all the systems we studied, the N = 400 data set is the least smooth and most difficult to equilibrate. To verify that we have obtained the true force−extension behavior for N = 400, we prepared three independent sets of simulations, all starting from different initial configurations. The results are shown as symbols of different colors in Figure 4d, which suggests that our force−extension results are reliable. In Figure 5, we plot together the force dependence of Rz from all chain lengths, normalized by N (Rz/N may be thought of as the average z component of the bond vectors). The data collapse to a universal force−extension law. This universal

behavior deviates from the Gaussian prediction (dashed line) at a force value of about 0.2 kT/σ. This deviation is not related to chain entanglements, as chains are topologically equilibrated in our simulations with the help of chain rebridging moves, but is caused by the fact that a physical chain has a finite length. To model this finite-extensibility effect, we have tried to fit the data with a freely joined chain model (FJC), which predicts that24 R z = Nl0 3(FlK /(kT ))

(21)

where 3 (x) ≡ coth(x) − 1/x is the Langevin function, l0 the mean bond length (l0 = 0.997σ in our model), and lK the Kuhn length defined by lK ≡ b2/l0. Because 3 (x) approaches x/3 for small x, eq 21 reduces to the Gaussian form at small force F. We recall the origin of the Kuhn length as the effective freely joined step size (nKlK = Nl0, and nKlK2 = Nb2, nK being the number of Kuhn steps). The Kuhn length lK = b2/l0 increases for stiff chains.

Figure 4. Extension-force relation for chains of N = 50 (a), 100 (b), 200 (c), and 400 (d). Dashed lines are predictions of linear response theory. For N = 400, data points of different colors are results from independent simulations. 1664


Macromolecules

Article

Evaluating eq 21 with the Kuhn length calculated by b2/l0 leads to the black curve in Figure 5. Compared to the Gaussian theory, eq 21 qualitatively captures the reduced extension at large force, but still overpredicts the simulated results. The simulation results behave as if the chains are more flexible, and hence less easily extended, once a force has been applied. We argue that the Kuhn length, which measures the length at which bond vectors are correlated, should be a function of F. In the undeformed state, lK is slightly greater than the bond length l0 (lK ≃ 2σ) because short distance repulsions prevent the (i + 1)st bead of a chain from back-folding onto the (i − 1)st bead. But in the strongly deformed state, since the chains have already been stretched, the beads are increasingly unlikely to make a backward move, so will feel increasingly smaller effects of repulsion from next neighboring beads. This implies that the Kuhn length in oriented melts should be lower than the unperturbed value. For chains close to the maximum extension, we expect the Kuhn length to be the same as the average bond length l0, since in that limit the chain configurations approach those of the freely joined chain. To incorporate this diminished effect of next-neighbor exclusions on lK for chains under tension, we propose the following form for the F dependence of the Kuhn length (lK,0 = 2.0σ and l0 = 0.997σ for our model): lK (F ) = lK ,0 + (l0 − lK ,0) tanh(αFl0/(kT ))

have been recorded over a total time of 62.5τe (τe is the Rouse time of an entanglement strand;1 the value was determined as reported in ref9.). Furthermore, throughout our work, the chain rebridging moves are used prior to the isoconfigurational ensemble simulation, so the system configurations even for stretched chains are well equilibrated, and the estimated tube diameter is that of an equilibrium ensemble. The tube diameter estimated from the isoconfigurational ensemble simulation shows a weak dependence on the duration τa of the MD trajectory, as discussed in our previous work.9,25 In Figure 6a, we plotted the dependence of tube diameter on

Figure 6. (a) Dependence of tube diameter on bead index in the undeformed state. Curves from bottom to top are obtained at τa = 40, 50, 100, 150, 200, and 250 frames, which corresponds to 10, 12.5, 25, 37.5, 50, and 62.5 τe respectively. (b) Dependence of characteristic ratios c1 ≡ (L/a)/(N/Ne) (left side) and c2 ≡ a/(Ne1/2b) (right side) on the averaging time τa.

(22)

The tanh function represents a crossover, such that (1) lK should decrease particularly for small forces and (2) lK should approach l0 for large forces. Here α is an adjustable parameter, obtained by fitting to our simulated results. The fitted result with α = 0.6 is shown as the blue solid curve in Figure 5, which represents well our data for the full range of force values. The total fractional variation (l0 − lK,0)/lK,0 in the Kuhn length is about 50%. With eq 22 and α = 0.6, lK(f) has shifted down by about 20% at our largest applied force F = 0.7kT/σ. Tube Diameter. We use the isoconfigurational ensemble simulation proposed in our previous work9,25 to estimate the tube diameter. We take a well equilibrated system configuration as the starting point, and run a sequence of MD simulations. All these MD simulations start from the same initial configuration (isoconfiguration), but each has a different set of bead velocities drawn from the Boltzmann distribution. The idea is to allow the beads to explore the tube confinement through different paths. In the present work, a key point in conducting isoconfigurational simulations is that the pulling force at chain ends continues to be applied. Otherwise, the isoconfigurational ensemble would describe chain retraction dynamics. For each MD trajectory in the isoconfigurational simulation, we record the bead positions periodically. The set of positions for a single bead recorded from all MD trajectories form a “cloud”. The shapes of the cloud show vividly how motions of the beads are constrained to lie near an underlying smooth curve, which we interpret as the tube primitive path. If we map each cloud onto a density function, the tube primitive path can be identified as the locus of the maximum density of the cloud.25 Once the tube primitive path has been identified, the extent of the transverse spreading around the path can be used as a measure of tube diameter. The details have been reported in refs 9 and 25. To improve statistics, for each force value, we conducted isoconfigurational simulations on six equilibrated system configurations. Each isoconfigurational ensemble contains 100 MD trajectories. For each trajectory, 250 evenly spaced frames

the bead index s in the undeformed state. The values are obtained by calculating the radius of gyration matrix of bead cloud with respect to the bead position on the primitive path, computing its eigenvalues, and averaging the square root of the two smaller eigenvalues (this is a radius; the diameter we use is twice this value). In addition, the values are averaged over the 10 molecules and six independent sets of simulations. Different curves in Figure 6a are results obtained at different τa (ranging from 40 frames (10 τe) to 250 frames (62.5 τe)). The lower curves (smaller tube diameter) correspond to smaller τa values. The origin of this dependence was addressed in refs 9 and 25. Here we focus on the bead index dependence. Figure 6a clearly shows that for beads in middle part of tube the tube diameters are nearly constant. For beads closer to the tube ends, the apparent tube diameter increases, which results from the chain end motions that can occasionally escape the tube confinement. As shown in the figure, this end effect persists for about 60 beads, close to the value of Ne. The apparent tube diameter for beads at the tube ends is higher than that in the middle by about 70%. In the following discussion, the values we used for the tube diameter are always the averages of values over the 200 beads in the middle of tube. In the standard tube theory,1 the tube primitive path contour is approximated as a freely jointed random walk of N/Ne steps of fixed length a. The primitive path contour length L is then given by L = (N/Ne)a. Each step consists of a strand of Ne beads executing a random walk, with a2 = Neb2. In our previous work on direct observations of the primitive path and tube using isoconfigurational averaging,9 we made several quantitative measures of the tube and primitive path properties, which we now reprise and analyze. The primitive path for a ring of 800 beads in a melt was found to have a contour length L of about 120σ. Furthermore, the primitive path was smooth, with a persistence length LP of about 40 monomers. Treating the primitive path itself as a semiflexible random walk, this implies a Kuhn segment for the primitive 1665


Macromolecules

Article

Figure 7. (a) Averaging time dependence of the tube diameter at different force values. From top to bottom, F = 0, 0.044, 0.091, 0.140, 0.190, 0.246, and 0.313 kT/σ, respectively. (b) Force dependence of the tube diameter at different averaging time τa. Circles: τa = 100 frames. Squares: τa = 250 frames. Closed symbols: a from gyration radii of monomer “clouds”; open symbols: from minimum distance to primitive path (see main text).

path NK equal to 2L P, or about 80 monomers. The corresponding Kuhn length can be computed as the meansquare end-to-end distance of this primitive path Kuhn segment, LK = (NKb2)1/2, equal to about 12.6σ. Alternatively, we can equate the tube contour length L to (N/NK)LK, again obtaining LK of about 12σ. We emphasize that N K can be interpreted as the entanglement length, in the sense that the primitive path scales as a freely jointed walk of N/NK steps of length LK, with LK2 = NKb2. Measures of the length scale for transverse fluctuations of the physical chain, such as the tube diameter a (twice the transverse radius of gyration of the clouds of monomer positions), will be proportional but not equal to LK. Now the physical chain is confined within the tube along the primitive path, and experiences a thermal tension FT sufficient to extend it to the proper length L. Writing the linear force− extension relation for the physical chain as R = (Nb2/(3kT))F and equating R to L = (N/NK)LK, we find the thermal tension FT as FT =

3kT LK

We see clearly the presence of two regimes: the tube diameter is nearly constant when F is less than 0.1 kT/σ, and decreases by up to about 25% for greater applied forces. The thermal tension computed from eq 23 is FT ≈ 0.25kT/σ, which means that the data of Figure 7b is clearly in the crossover region. To describe these results quantitatively therefore requires not simply a scaling argument, but a theory of the crossover from weak to strong chain tension, which we shall provide in the next subsection. We have also extracted values for tube diameter by a slightly different method, designed to mitigate the apparent timedependence of the tube diameter and primitive path length that results from clouds of monomer positions tending to curve around sharp corners of the primitive path. This “cornercutting” leads to the center of the cloud being located “inside the bend” in the primitive path. In the alternative method, we identify the primitive path as the locus of the points of maximum density in the cloud of positions for each monomer. Then, we define the distance from each point in the cloud to the “nearest point” on the primitive path. The search for the nearest point on the path is local, in that we start with the point corresponding to the cloud maximum density, and search for the nearest local minimum in the distance from the cloud point to the path. Once this minimum distance is found, it is treated as a datum point for the histogram of the tube radius. The tube radius is then computed as an average over points in each cloud, monomers along the chain, and configurations among the isoconfigurational average, for each value of the applied tension. The results are shown in Figure 7 (open symbols) for two different sampling times, τa = 100 (circles) and τ = 250 (squares). Unfortunately, this heuristic approach does not eliminate the modest time-dependence of the apparent tube diameter. Crossover of Tube Diameter. The key step in developing the various forms of the entanglement ansatz is estimating the dimensions of the volume pervaded by a strand of Ne beads. In the scaling arguments, the pervaded volume V is estimated differently for isotropic melts and chains under tension. For isotropic melts, it is estimated as a3; for chains under tension, it is estimated using the volume of the cylinder pervaded by a test strand, as described in the section Entanglement Ansätze for Chains under Tension. To proceed beyond separate scaling arguments for these two regimes, we need a crossover function for the pervaded volume V, which we propose to be

(23)

Note that LK ≈ 12σ is about twice the tube diameter a ≈ 6σ, itself defined as twice the gyration radius of the cloud of monomers surrounding the primitive path. This factor of 2 is a consequence of the semiflexibility of the tube, and is important for quantitative determination of the thermal tension. We emphasize that different detailed definitions of the entanglement length and tube diameter will give different numerical values, with factors of order unity between different definitions. For example, the chain shrinking algorithm gives Ne = 67;9 as described above, the primitive path persistence length gives Ne = 80;9 and the topological entropy counting gives a value between 52 and 67.16 We can as well obtain a value for Ne from the original Lin− Noolandi relation for isotropic melts, in the form eq 3; with b2 ≈ 2σ2 and ρ = 0.7/σ3 so that Ω = σ3/0.7, we find Ne = 128. A useful discussion of some aspects of these ubiquitous factors of order unity has been given in ref 26. In Figure 7a, we plotted the τa dependence of tube diameter for different forces. Each curve corresponds to one force value. Despite the dependence on τa, it is still apparent that the tube diameter decreases with force. To minimize the dependence on the averaging time, in Figure 7b, we plotted the force dependence of the tube diameter normalized by the undeformed values. The difference between results at two different τa values is indeed small.

V ≡ (det .)1/2 1666

(24)


Macromolecules

Article

30. In the comparison, we first adopt the value for Ne0 of NK = 80 monomers, and then regard this value as adjustable, to see whether a reasonable value can account for the simulation results. The curves in Figure 8 are predictions of eq 30 under ansatz I, with Ne0 = 80 (solid) and Ne0 = 60 (dashed) respectively, corresponding to values for FT of 0.25kT/σ and 0.29kT/σ respectively.

where . is the radius of gyration tensor of a strand of N beads. The components of . are defined by Gαβ ≡

1 N

N

∑ ⟨(R i − R c)α (R i − R c)β ⟩

(25)

i

where Ri is position vector of the i-th bead, Rc ≡ is the center of mass of the strand, and α and β are Cartesian indices. For strands obeying Gaussian statistics, all off-diagonal terms in . vanish. The diagonal terms are found to be (1/N)∑Ni Ri

Gxx = Gyy = Gzz =

Nb2 18

2 Nb2 N 2b4 ⎛ F ⎞ Nb2 ⎛⎜ 3 ⎞ ⎜ ⎟ = + 1 + nf 2 ⎟ ⎝ 18 12 ⎝ 3kT ⎠ 18 2 ⎠

(26)

The Gxx and Gyy components come from fluctuations of the bead-to-bead separation vectors. The Gzz component includes contributions from both fluctuations and the average drift of successive monomers along the z direction as a result of the applied tension. In the second line, we have normalized F by the thermal tension FT = 3kT/a; the ratio N/Ne0 is denoted by n. As we have different numerical estimates of Ne0, the corresponding value of FT and hence the detailed behavior of the crossover will depend on which estimate we adopt. For example, we may settle on NK as described above (NK ≈ 80 for our bead−spring chain) as the numerical value of Ne0, whereupon FT = 3kT/LK ≈ 0.25kT/σ is the corresponding value of thermal tension. We shall discuss this choice further below. Ansatz I (Lin−Noolandi: fixed number of strands in the pervaded volume of an entanglement strand) can be written as V (Ne ,0 , 0) V (Ne , F ) = = C′ Ne Ω Ne ,0 Ω

Figure 8. Symbols: tube diameter simulation results (same as Figure 7). Curves are results of ansatz I theory, with Ne0 = 80 (solid) and Ne0 = 60 (dashed).

Next we consider ansatz II (Colby−Rubinstein−Milner: fixed number of contacts in the pervaded volume of an entanglement strand). This can be written (see eq 4) V (Ne , F ) p3

(27)

ne3(1 + (3/2)nef 2 ) = 1

(28)

(1 + 6f 2 )1/2 − 1 3f 2

(29)

For small f, ne approaches unity; for large f, we find ne approaches (2/3)1/2/f, in agreement with our scaling argument for ansatz I. The tube diameter follows from Gxx or Gyy evaluated at Ne = Ne,0ne, which measures the transverse fluctuation. So we have a = (18Gxx)1/2 = ne1/2a0

(32)

with the relation between a(F) and ne(F) given as before by eq 30. Ansatz II evidently gives weaker dependence for Ne(F) for large tension, with Ne(F) going as F−1/2, consistent with our earlier scaling argument. Recall that ansatz II gives reasonable results for diluted melts, although more recent comparison to rheological data as well as our own simulation results for entanglement of diluted melts is more consistent with Ne(ϕ) scaling as ϕ−1 than with ϕ−4/3 as predicted by ansatz II. As for ansatz I, we compare the predicted force dependence of the tube diameter ratio a(F)/a0, first adopting the value Ne0 of NK = 80 monomers, then adjusting this value for a better fit. The curves in Figure 9 are predictions of eq 30 under ansatz II, with Ne0 = 80 (solid) and Ne0 = 124 (dashed) respectively, corresponding to values for FT of 0.25kT/σ and 0.20kT/σ respectively. The weak dependence under ansatz II of a(F) on tension requires a considerably larger value of Ne (hence relatively more sensitivity to small applied tension) to be consistent with the simulation results. Finally, we develop a crossover formula for ansatz III (present work: fixed topological complexitymore specifically, fixed number of contacts per entanglement strand, when deformed to an isotropic melt at fixed topological complexity). The scaling assumptions for this ansatz amount to (see eq 15)

The physical solution of this quadratic is ne =

(31)

Again enforcing that Ne,0 obtained here is the same as for isotropic melts fixes C′ according to C′3 = 18−3/2C3. Using eq 26, we find that ne = Ne(F)/Ne0 now satisfies

Here for clarity we have shown the dependence of V on Ne and F explicitly. Enforcing that Ne,0 obtained in this way equals the result for isotropic melts from eq 1 fixes C′ to be C′ = 18−3/2C. Using eq 26, we find that ne = Ne(F)/Ne0 satisfies ne(1 + (3/2)nef 2 ) = 1

= C′3

(30)

where the factor 18 is added for consistency with eq 2, and a0 is the unperturbed tube diameter. We have argued that ansatz I does not give reasonable results for diluted melts, and so is a doubtful choice for describing entanglement in melts under tension. For completeness, we compare its predictions for melts under tension to our simulation results for the force dependence of the tube diameter, expressed as the ratio a(F)/a0, using eqs 29 and 1667


Macromolecules

Article

Marrucci and de Cindio27,28 proposed that the tube diameter scale with F−1/2, which was also used by Wagner et al.29,30 They argued by assuming a “constant tube segment volume”,27 a(l′)2l′ = a(l)2l, which implies that the tube diameter scales with the ratio (l′/l)−1/2, where l′ and l are lengths of stretched and unstretched tube segments. The ratio l′/l measures the extent of stretch. Theories built upon this assumption have accounted for some aspects of polymer melt rheology under extensional flow sufficiently fast to stretch chains (e.g., refs.29,30 and28). The GLaMM model8 likewise makes an assumption regarding the tube for chains under tension. Ref 8 considered two possibilities: (1) the number of entanglement points along a given chain remains fixed under affine deformation; or (2) the primitive path length between entanglement points remains fixed. GLaMM makes the latter assumption, under which the number of entanglement points increases as the chain stretches (Ne decreases). The predictions of the GLaMM model for nonlinear rheology agree rather well with a wide range of experimental data. The scaling found under ansatz I (Ne( f) scaling as f−1) is the same as was assumed in the GLaMM model of nonlinear rheology of entangled melts. Unfortunately, ansatz I does not give reasonable results when applied to diluted melts (Ne(ϕ) scaling as ϕ−2). Ansatz II gives reasonable results for diluted melts (Ne(ϕ) scaling as ϕ−4/3), but too weak a dependence of Ne(f) on applied tension (Ne( f) scaling as f−1/2). Ansatz III gives both a reasonable result for diluted melts (Ne(ϕ) scaling as 1/ϕ), and a reasonable result for melts under tension (Ne(F) scaling as F−2/3). While this is not as strong as assumed in the GLaMM model, it seems possible that this scaling could (with the adjustable parameters necessary in GLaMM) be made to give reasonable results for nonlinear rheology. Qualitatively, the GLaMM assumption fights against the possibility of a stress maximum, since the number of entanglement segments increases as the chains become stretched. This means that a factor of (R’(s)/a)2 times many tube hops are required to relax a chain segment. Also, the spring constant of each entanglement strand is increased by a factor of |R’(s)|/a. The net effect is a convective constraint release (CCR) term of the form8

Figure 9. Symbols: tube diameter simulation results (same as Figure 7). Curves are results of ansatz II theory, with Ne0 = 80 (solid) and Ne0 = 124 (dashed).

λ 2 = Gzz /Gxx N′p = Np/λ Ne/N′p = C 2

(33)

Combining these results, we find ne = 1/λ; using eq 26, we find ne2(1 + (3/2)nef 2 ) = 1

(34)

with a(F) again given in terms of ne(F) by eq 30. Ansatz III gives dependence for Ne(F) for large tension intermediate between ansatz I and III, with Ne(F) going as F−2/3. We again compare the predicted force dependence of the tube diameter ratio a(F)/a0, first with Ne0 = NK = 80 monomers, then adjusting for a better fit. The curves in Figure 10 are predictions of eq 30 under ansatz III, with Ne0 = 80

∂R(s) 3ν a R″(s) + g (̃ s , t ) = ... + 2 |R′(s)| ∂t

Figure 10. Symbols: tube diameter simulation results (same as Figure 7). Curves are results of ansatz III theory, with Ne0 = 80 (solid) and Ne0 = 88 (dashed).

(35)

The above equation describes the time evolution of the tube contour R(s): s parametrizes the segments along the tube contour, ν is the constraint release rate, and g̃ is the random noise needed by the thermal equilibration. If the reduction in Ne for chains under tension were less pronounced than assumed in GLaMM, it may be that a stress maximum would reappear in the predictions of the theory. Effects of Nonlinear Extension. Our simulation results are limited to rather modest forces, because of the difficulty of equilibrating stretched systems at large tension. But for our bead−spring chains, the Kuhn length of the physical chain lK is about 2σ, which means the characteristic force for nonlinear extension is FK = kT/lK = 0.5kT/σ. This is not much larger than the force FT to perturb the chain within an entanglement strand, estimated above as about 0.25kT/σ. In other words, the Kuhn length of the physical chain for our bead−spring model is not much smaller than the tube diameter.

(solid) and Ne0 = 88 (dashed) respectively, corresponding to values for FT of 0.25kT/σ and 0.24kT/σ respectively. The crossover function of ansatz III accounts for the simulation results with only very modest adjustment of the value of Ne0 away from the value we extract from treating the primitive path for our bead−spring chains as an effective freely jointed path of step size given by the Kuhn segment for the primitive path. Relation to Previous Work. The scaling of the tube diameter and entanglement molecular weight with tension is a key element of contemporary models of nonlinear rheology. The F−1/3 dependence of the tube diameter we obtained from version III of the Lin−Noolandi ansatz differs quantitatively from the assumptions made by previous works concerned with nonlinear rheology. 1668


Macromolecules

Article

nef2 (x)2 (f1 (x) + (3/2)nef 2 (f2 (x))2 ) = 1 (ansatz I)

This is typically the case for chains that are relatively stiff and not very bulky; real polymers of this sort include polyethylene (PE) and poly(ethylene oxide) (PEO). It is therefore of interest to ask what happens to the tube diameter at sufficiently large applied tension that the chains are close to maximum extension. For large enough tension, even the packing blobs are perturbed. However, the packing length p is typically much smaller than the Kuhn length, so there is a considerable range of tensions for which we may regard the packing blobs as unperturbed, which simplifies our analysis. We recalculate the gyration radii Gxx = Gyy and Gzz for a chain under strong tension by modeling the chain as a freely jointed chain with step size equal to its Kuhn length lK. Here, we ignore the complication that the Kuhn length for our bead− spring chains appears to shift to smaller values at sufficiently large applied force; lK decreased only by about 20% at the largest force we applied (0.7kT/σ, in excess of FK = 0.5kT/σ). The gyration radii for a strand of N beads using the FJC model can be calculated as Gxx = Gyy =

Gzz =

ne3f2 (x)2 (f1 (x) + (3/2)nef 2 (f2 (x))2 ) = 1 (ansatz II) ne2(f1 (x) + (3/2)nef 2 (f2 (x))2 )/f2 (x) = 1 (ansatz III) (40)

These results evidently reduce to those of the previous section for small x (applied force less than the Kuhn force). By using eq 30, we find that regardless of which ansatz is assumed, the tube diameter satisfies the relation ⎛ 18G ⎞1/2 a xx ⎟ = ⎜⎜ = (nef2 (x))1/2 2⎟ a0 N b ⎝ e ,0 ⎠

For the range of forces for which we have simulation data, the above crossover formulas result in relatively small corrections to the corresponding formulas without nonlinear elasticity effects, certainly in comparison to the changes that result from adjusting the value of Ne0. However, for applied tensions much larger than the thermal tension, for which f1(x) and f 2(x) approach their large-x limits, the scaling of Ne(F) is strongly altered by the nonlinear effects. For either ansatz I or II, Ne(F) exhibits a minimum, at an applied tension of about 1.3 times the Kuhn force (i.e., F ≈ kT/ σ). For large force, Ne(F) grows as F for ansatz I and as F1/2 for ansatz II. Thus, if either of these two assumptions were physically reasonable, chains near full extension would be very weakly entangled. For chains of finite length N, the system would eventually disentangle when Ne grew to a value comparable to N. Under either of these two assumptions, the shape of Ne/Ne,0 as a function of pulling force can be understood qualitatively as follows. For small tensions, chains freely wander into the volume occupied by neighboring chains. Ne decreases for small tension because the pervaded volume dimension increases along the pulling direction. For large tensions, the extent of transverse fluctuations is increasingly limited, and the strand length needed to produce an entanglement increases accordingly. For ansatz III, Ne(F) does not exhibit a minimum, but does weaken its dependence on F for large tensions, to F−1/3. If we were able to obtain reliable results for the tube diameter for applied forces well in excess of the Kuhn force (but less than kT/p, so that packing blobs are unperturbed), we could presumably distinguish the presence or absence of a minimum entanglement length. Because of the nonlinear elasticity effects, the corresponding scalings for the tube diameter no longer agree with the Gaussian behavior in which a scales as Ne1/2. For large tensions, we find scaling results for a of F0 (ansatz I), F−1/4 (ansatz II), and F−2/3 (ansatz III). The tube narrows under the assumptions of each ansatz (reaching a finite limit (ansatz I) or continuing to decrease (ansatz II and III)) whether or not Ne displays a minimum, because of the strong suppression of transverse fluctuations in the stretched chains. Table 1 summarizes the scaling predictions of ansatz I, II, and III for the scaling of the entanglement molecular weight Ne with volume fraction ϕ and applied tension F. For most polymers, the Kuhn length is several times larger than the packing length (see Table 2). Thus, the force kT/lK required to bring most polymers to near full extension is

(N /NK )lK 2 ⟨sin 2 θ ⟩, 12

(N /NK )lK 2 (N /NK )2 lK 2 ⟨(Δcos θ )2 ⟩ + ⟨cos θ ⟩2 6 12 (36)

where we have used the definition ⟨(Δ cos θ) ⟩ ≡ ⟨cos θ⟩ − ⟨cos θ⟩2. (Here NK is now the number of monomers in a physical Kuhn segment, of length lK; a strand of N monomers consists of N/NK Kuhn segments.) The averages are evaluated using the probability weight P(cos θ) ∝ exp(FlK cos θ/(kT)) for −1 ≤ cos θ ≤ 1, which yields 2

⟨cos θ ⟩ = 3(x)⟨sin 2 θ ⟩ = (2/x)3(x)

2

(37)

with 3 (x) the Langevin function given by 3 (x) = coth(x) − 1/x, and x is the dimensionless force F/FK, with FK the “Kuhn force ” given by FK = kT/lK. Applying these results to an entanglement strand of length N = Ne, we can write Gxx = Gyy = Gzz =

Neb2 f (x ) 18 2

Neb2 (f (x) + (3/2)nef 2 (f2 (x))2 ) 18 1

(38)

in which we have defined f1 (x) = 3(1 − (2/x)3(x) − 32(x)) f2 (x) = (3/x)3(x)

(41)

(39)

Both f1(x) and f 2(x) approach unity for small x; for large x (large applied tension), f 2(x) approaches 3/x with exponentially small corrections while f1(x) approaches 3/x2. Making use of these results for Gxx and Gzz, the crossover equations for ansatz I, II, and III take the form 1669


Macromolecules

Article

Table 1. Scaling of ne = Ne/Ne,0 with Volume Fraction ϕ and Tension F, under Ansatz I, II, and IIIa ansatz

solution

F < FT

F > FT

F ≫ FT

I II III

ϕ−2 ϕ−4/3 ϕ−1

1 1 1

F−1 F−1/2 F−2/3

F F1/2 F−1/3

smaller than a, but will begin to be comparable to p (e.g., for PDMS a/lK ≈ 6, but lK/p ≈ 3). So polymers will typically either display some effect of semiflexibility of chain segments within their tube (small a/lK),31 or some force-dependence of the Kuhn length (small lK/p).

■

SUMMARY We have studied the effect of chain orientation and stretching on the tube diameter for polymer melts, both theoretically and by simulation methods. The chain orientation is introduced by applying pulling forces of opposite directions at the chain ends. In the absence of pulling force, chain configurations are well described by the random walk statistics down to the scale of the packing length p = Ω/b2, which as implied by the Lin− Noolandi relation is about five percent of the tube diameter.7 In the presence of pulling force, the chain configuration can be thought of as a string of tension blobs, the blob size of order ξ = 3kT/F. In the small force regime, the tension blob size is larger than the tube diameter, so that the configurations of entanglement strands are unperturbed. Once the tension blob size falls below the tube diameter, the force perturbs the configuration of the original entanglement strand, whose pervaded volume can then be considered as a cylinder. It is not obvious how to develop an entanglement ansatz to describe chains in a melt under tension. For isotropic melts, there are three equivalent formulations: (1) fixed number of neighboring strands cohabiting the volume V pervaded by an entanglement strand (Lin−Noolandi, ansatz I); (2) fixed number of close contacts between chains (of volume p3) in the volume V (Colby−Rubinstein−Milner, ansatz II); (3) fixed number of close contacts between an entanglement strand and neighboring chains (present work, ansatz III). These three formulations are not equivalent for diluted melts (such as mixtures of long entangled chains with shorter, quickly relaxing chains) or for melts of chains under tension. For diluted melts, Ne(ϕ) scales as ϕ−2 (ansatz I), ϕ−4/3 (ansatz II),7 or ϕ−1 (ansatz III). Ansatz I is not consistent with rheological data on the concentration dependence of plateau modulus. Ansatz II was believed to be consistent with such data,7 but more recent studies suggest that ansatz III better accounts for the observed behavior.11 We have extended all three formulations to the case of chains in a melt under tension. The extension of ansatz III is subtle, and may be briefly described as “fixed number of contacts per entanglement strand, when deformed to an isotropic melt at fixed topological complexity” (see subsection Ansatz III above). We find Ne(F) scales as F−1 (ansatz I), F−1/2 (ansatz II), and F−2/3 (ansatz III). The tube diameter a(F) is predicted to scale as Ne1/2 so long as the applied force is not large enough to bring the chains to near full extension. To test these predictions, we simulated a melt of chains oriented by pulling the chains at the ends, and equilibrated topologically by various Monte Carlo moves that permit chains to cross. We first studied the force−extension law in these systems, and found the system response becomes nonlinear when the chain is extended by only about 10% of its fully extended length. The resulting force−extension relation can be described by a modified Langevin function, for which the Kuhn length lK depends weakly on the external force. We then performed isoconfigurational ensemble simulations for the stretched chains, and estimated the tube diameter by the width of the monomer fluctuations transverse to the primitive

a

In solution and for small tension, a/a0 scales as ne1/2; for large tension, a/a0 scales as (ne f 2(x))1/2 [eq 41].

Table 2. Mean-Square End-to-End Distance per Mass R2/M; Density, ρ; Tube Diameter aexp;3 Kuhn Length lK;24 and Predicted Tube Diameter, acalc [Eq 2 ] for Common Polymers at 140 °Ca polymer

R2/M

ρ

M0

lK

p

aexp

acalc

PE PEO PMMA PS PP PB PDMS PI

1.25 0.805 0.425 0.434 0.670 0.876 0.457 0.625

0.784 1.06 1.13 0.969 0.791 0.826 0.895 0.83

14.0 44.1 100.1 104.2 42.1 54.1 74.2 68.1

14 11 17 18 11 9.6 13 8.2

1.69 1.94 3.46 3.95 3.13 2.29 4.06 3.20

33 38 67 77 61 44 79 62

36 40 75 85 51 70 94 75

a Packing lengths p from M0/(b2ρ). Key: PE, polyethylene; PEO, poly(ethylene oxide); PMMA, poly(methyl methacrylate); PS, polystyrene; PP, poly(propylene); PB, 1,4-polybutadiene; PDMS, poly(dimethylsiloxane); PI, 1,4-polyisoprene. Lengths in Å; masses in g/mol.

significantly smaller than the force required to perturb chain on the scale of the packing length lp. Roughly speaking, chain stiffness can arise from the intrinsic stiffness of the backbone (from relatively large trans−gauche energy differences, or small deflection angles between successive dihedral bonds), or it can arise from repulsive interactions between successive monomers and their side groups along the chain. For unusually bulky polymers, it may happen that the Kuhn length becomes more comparable to the packing length. For such chains, the Kuhn length may be determined in part by repulsions between successive monomers. Then, a force large enough to perturb the chain on the scale of the Kuhn length may also begin to reduce intrachain monomer−monomer repulsion, thus reducing the effective Kuhn length. We observed a modest dependence of the Kuhn length on applied force in the elastic response for our bead−spring model (see eq 22 and surrounding discussion). For our bead−spring model, the packing length p = Ω/b2 is about 0.7σ, while the Kuhn length lK is about 2σ; so lK/p is about 2.8. The more bulky, flexible chains in 2 such as PB, PDMS, and PI have ratios of lK/p close to 2.7. We may therefore expect these polymers to exhibit anomalous stretching behavior similar to our simulation model. The other polymers listed in Table 2 all have much higher values of lK/p, so their Kuhn lengths are expected to be independent of force. By definition, for flexible entangled polymers the Kuhn length lK must lie somewhere between p and a. Since the ratio a0/p is fixed according to the Lin−Noolandi relation, there is a fixed range of length scales in which lK must reside. If a polymer (such as PE or PEO) is rather stiff and skinny, lK will begin to approach a (a/lK for PE is about 2.5). Or, if a polymer (such as PB, PDMS, or PI) is rather bulky and flexible, lK will be rather 1670


Macromolecules

Article

so it seems reasonable to argue that the effect of tension on the tube in both cases is similar. Although there might be some differences, we expect the stretched system to be a good proxy for systems under flow. The deformation dependence of the tube diameter in topologically quenched systems such as polymer gels is quite different; we address this system in a separate paper.25

path. Because of difficulty in equilibrating strongly stretched chains, our results are limited to applied forces of order the thermal tension, so that the entanglement strands are only modestly stretched. Thus, our results for a(F) are in a crossover region, and do not display simple scaling dependence on force. To distinguish the behaviors of ansatz I, ansatz II, and ansatz III with regard to a(F), we compute for each formulation a crossover function. We determine the necessary parameters by analyzing the dimensions of conformations of isotropic chains and their primitive paths. Any of the three formulations can be made to account for the simulation results if we adjust the isotropic melt value of Ne for our bead−spring model. However, ansatz I and ansatz II require substantial adjustments of Ne away from the value extracted from the Kuhn length of the semiflexible primitive path. Only ansatz III is consistent with the data without large adjustments of Ne. Since only ansatz II and ansatz III are consistent with the concentration dependence of Ne(ϕ) for diluted melts, this points to ansatz III as the proper formulation of the original Lin−Noolandi idea. While ansatz IIfixed number of close contacts within the pervaded volume of an entanglement strandsuggests a “collective” origin for entanglements, ansatz IIIessentially, a fixed number of close contacts per entanglement strandpoints toward a more “binary” origin of entanglements. We note that the scaling dependence of Ne(F) predicted by ansatz III is somewhat weaker than that assumed in the successful GLaMM model of the nonlinear rheology of linear chains. However, the GLaMM model contains an important adjustable parameter (which controls the strength of constraint release effects); it may be that a version of the GLaMM model with weaker scaling for Ne(F) may also perform well with respect to experimental data. In principle, the force dependence of tube diameter could also be investigated with the chain-shrinking algorithms.32−35 But we suspect that the moderate fractional change in tube diameter (25%) is too small to be reliably identified with those methods, in which the entanglement length is revealed by counting binary contacts of the shrunken chains. Since the contacts are relatively few in number, measuring a small change in Ne would be statistically challenging. Our simulations are limited to relatively small force values. For applied forces above the “Kuhn force” (kT per Kuhn length), chains begin to exhibit nonlinear elastic response, and approach full extension. We have extended our crossover functions to describe this regime as well, in which the scaling of Ne(F) changes because of the nonlinear elastic effects. For ansatz I and ansatz II, Ne(F) is predicted to exhibit a minimum, at a force of about 1.3 times the Kuhn force; for ansatz III, Ne(F) continues to decrease, albeit more weakly as F−1/3. For all three formulations, the tube diameter decreases monotonically with increasing force. Our model is limited to forces less than kT/p, such that the packing “blobs” are unperturbed. If we could manage to obtain simulation results for larger applied forces, it may be possible to distinguish ansatz II from ansatz III unambiguously by the presence or absence of a minimum in Ne(F). The system we studied here is oriented by pulling chain ends, which differs from the system oriented by applying a shear flow but has the advantage of being generated from a well-defined and topologically equilibrated ensemble. In systems under flow, chains certainly cannot cross, but over time they do explore different topological configurations, and they are under tension,

■

APPENDIX Here we briefly reprise the scaling results of ref 7 for dependence of the entanglement length Ne(ϕ), under the assumptions of ansatz I (constant number of neighboring strands in the pervaded volume V of an entanglement strand) and ansatz II (constant number of binary contacts in the pervaded volume V). For ansatz I, following eq 1 we write ϕa3(ϕ) =C Ne(ϕ)Ω

(42)

in which the factor of ϕ appears because only a fraction ϕ of the pervaded volume actually contains neighboring strands. We use the relation a2 = Neb2 to eliminate a(ϕ) in favor of Ne(ϕ), and then divide the above equation by its ϕ = 1 limit to obtain Ne(ϕ) = ϕ−2 Ne0

(43)

For diluted systems, ansatz I amounts to the assumption that each polymer strand is covered with a uniform “coating” of diluent, i.e., that the monomer volume Ω has been replaced by Ω/ϕ. Actually, a Θ solvent is distributed more or less randomly in a diluted melt, which is why ansatz I is not reasonable to describe diluted melts. In ref 7, ansatz II was developed to describe diluted melts. The pervaded volume of an entanglement strand was assumed to contain a fixed number of binary contacts between “packing blobs” of volume p3 (see eq 4). In a diluted melt this takes the form ϕ2a3(ϕ) = C3 p3

(44)

The factor ϕ2 arises by breaking up the pervaded volume into sites of volume p3, only a fraction ϕ of which are occupied, at random; the number of contacts is then reduced by a factor ϕ2. Once again replacing a(ϕ) in terms of Ne(ϕ), and dividing the above equation by its ϕ = 1 limit, we obtain

■

Ne(ϕ) = ϕ−4/3 Ne0

(45)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Daniel Read for a close and patient reading of this paper in multiple versions, and many helpful comments and suggestions. We acknowledge support from the NSF (DMR1671


Macromolecules

Article

0907370) and the Donors of the Petroleum Research Fund, administered by the American Chemical Society (ACS-PRF 49964-ND7).

■

REFERENCES

(1) Doi, M.; Edwards, S. F. The Theory of Polymer Dynamics; Oxford University Press: Oxford, U.K., 1986. (2) McLeish, T. C. B. Adv. Phys. 2002, 51, 1379−1527. (3) Fetters, L. J.; Lohse, D. J.; Richter, D.; Witten, T. A.; Zirkel, A. Macromolecules 1994, 27, 996−998. (4) Lin, Y.-H. Macromolecules 1987, 20, 3080−3083. (5) Kavassalis, T. A.; Noolandi, J. Phys. Rev. Lett. 1987, 59, 2674− 2677. (6) Colby, R. H.; Rubinstein, M. Macromolecules 1990, 23, 2753− 2757. (7) Milner, S. T. Macromolecules 2005, 38, 4929−4939. (8) Graham, R. S.; Likhtman, A. E.; McLeish, T. C. B.; Milner, S. T. J. Rheol. 2003, 47, 1171−1200. (9) Bisbee, W.; Qin, J.; Milner, S. T. Macromolecules 2011, 44, 8792− 8980. (10) Adam, M.; Delsanti, M. J. Phys. (Paris) 1984, 45, 1513−1521. (11) Wang, Z.; Chen, X.; Larson, R. G. J. Rheol. 2010, 54, 223−260. (12) Manuscript in preparation. (13) de Gennes, P. G. Scaling concepts in polymer physics; Cornell University Press: Ithaca, NY, 1979. (14) Kremer, K.; Grest, G. S. J. Chem. Phys. 1990, 92, 5057−5086. (15) Weeks, J. D.; Chandler, D.; Andersen, H. C. J. Chem. Phys. 1971, 54, 5237−5247. (16) Qin, J.; Milner, S. T. Soft Matter 2011, 7, 10676−10693. (17) Fried, H.; Binder, K. J. Chem. Phys. 1991, 94, 8349−8366. (18) The software package is hosted at http://dmorse.github.com/ simpatico/index.html. (19) Morse, D. C.; Chung, J. K. J. Chem. Phys. 2009, 130, 224901. (20) Perkins, T. T.; Smith, D. E.; Larson, R. G.; Chu, S. Science 1995, 268, 83−87. (21) Li, L.; Larson, R. G. Macromolecules 2000, 33, 1411−1415. (22) Dobrynin, A. V.; Carrillo, J.-M. Y.; Rubinstein, M. Macromolecules 2010, 43, 9181−9190. (23) Smith, J. S.; Bedrov, D.; Smith, G. D.; Kober, E. M. Macromolecules 2005, 38, 8101−8107. (24) Rubinstein, M.; Colby, R. Polymer Physics; Oxford University Press: Oxford, U.K., 2003. (25) Qin, J.; So, J.; Milner, S. T. Macromolecules 2012, 45, 9816− 9822. (26) Larson, R. G.; Sridhar, T.; Leal, L. G.; McKinley, G. H.; Likhtman, A. E.; McLeish, T. C. B. J. Rheol. 2003, 47, 809−818. (27) Marrucci, G.; de Cindio, B. Rheol. Acta 1980, 19, 68−75. (28) Marrucci, G.; Ianniruberto, G. Macromolecules 2004, 37, 3934− 3942. (29) Wagner, M. H.; Schaeffer, J. J. Rheol. 1992, 36, 1−26. (30) Rolón-Garrido, V. H.; Wagner, M. H.; Luap, C.; Schweizer, T. J. Rheol. 2006, 50, 327−340. (31) Levine, A. J.; Milner, S. T. Macromolecules 1998, 31, 8623−8637. (32) Everaers, R.; Sukumaran, S. K.; Grest, G. S.; Svaneborg, C.; Sivasubramanian, A.; Kremer, K. Science 2004, 303, 823−826. (33) Zhou, Q.; Larson, R. G. Macromolecules 2005, 38, 5761−5765. (34) Kröger, M. Comput. Phys. Commun. 2005, 168, 209−232. (35) Tzoumanekas, C.; Theodorou, D. N. Macromolecules 2006, 39, 4592−4604.

1672


Tube Diameter of Oriented and Stretched Polymer Melts

Recommend Documents