A multichain polymer slip-spring model with fluctuating number of entanglements for linear and nonlinear rheology

Polymeric materials exhibit a wide spectrum of relaxation times that can range from picoseconds to milliseconds or even minutes.^1,2 As such, coarse-grained (CG) approaches are necessary to describe their long-time dynamics.³ It is well known that coarsening leads to a reduction of the effective interactions that are used to describe the physics of model, many-body systems. Interaction potentials become soft, i.e., the interaction energy between two interacting sites remains finite in the limit r → 0, where r is the distance between two CG particles. As a consequence, CG particles can cross each other, and important effects can be lost. This is particularly critical for polymeric materials, where the local non-crossability constraints that give rise to the peculiar dynamics of polymeric liquids can simply be lost. Indeed, the concept of “entanglements” plays a central role in the polymer physics literature.^4–7 Entanglements, which restrict the motion of high-molecular-weight chains, are introduced to explain the scaling of diffusion or viscosity with molecular weight. Highly coarse-grained models that rely on soft potentials lack this feature and are therefore unable to describe the rheology of high molecular weight materials.

Recent attempts to address this problem, including several by our own group, have adopted coarse-grained representations in which explicit entanglements are introduced at the two-molecule level.^8–13 In a many-chain representation of a polymer, the topological effects associated with non-crossability are described by effective fluctuating interactions, mediated by slip-springs (SS) that connect different macromolecules. The slip-springs are conceived in a manner that does not alter the structure and thermodynamics of the material at equilibrium, but they have a profound influence on dynamics. In our own past implementations,^11,12 which we refer to as “Theoretically Informed Entanglement Polymer Simulations,” or TIEPOS, it was shown that it is possible to reproduce key linear-response features that are characteristic of entangled polymers, including a scaling relation between relaxation time and molecular weight with an exponent larger than 3 and a plateau modulus. Such implementations were also shown to reproduce a number of experimentally observed non-linear rheology effects, in both steady shear and during start-up of shear flow.¹² Importantly, past work indicates that it is possible to achieve a quantitative description of experimental data for the linear rheology of monodisperse and bidisperse polystyrene blends.¹¹ 

On the other hand, at even coarser levels of description, single-chain formalisms have also been proposed to describe the rheology of entangled polymers.^{14–18,37–39} Among the latter models, the so-called discrete slip-link model (DSM)^19–22 has been shown to be capable of quantitative predictions for both linear and nonlinear rheologies.

As encouraging as past implementations of slip-springs in many-chain systems have been, a majority of them have implicitly assumed that the number of entanglements in the system is constant (one exception being Chappa et al.⁹) and does not change with applied fields (e.g., flow). Second, our own version of TIEPOS assumed that slip-springs diffuse in a continuous manner, according to dynamic equations and the corresponding diffusion coefficient. The “continuous” movement along polymer backbone requires careful consideration of the thermodynamic issues regarding the contribution of virtual springs. One key feature of the DSM is the use of a chemical potential bath for controlling the entanglement number, thus allowing this quantity to fluctuate and deviate from equilibrium values when non-equilibrium conditions are imposed. Furthermore, in DSM, slip-links perform discrete jumps along polymer chains and are forbidden to cross each other. Building on that body of work, we revisit the assumptions in our previous approach and present an implementation of TIEPOS where the number of entanglements is allowed to fluctuate, and where their motion is dictated by a jump process. For homogeneous materials, we find that the proposed TIEPOS approach is consistent with the predictions of single-chain theories, thereby raising intriguing prospects for systematic coarse-graining of polymer models for prediction of rheology in structurally homogeneous and heterogeneous materials such as composites, blends, or copolymers. We also show that TIEPOS are in quantitative agreement with experimental rheology data for polyisoprene melts.

For the sake of generality, we present the model in the context of an A-B polymer blend. The melt consists of N_c = n_A + n_B chain molecules in a volume V at temperature T. The macromolecules are represented by flexible linear chains described by the discretized Gaussian chain model. Each chain is composed of N_α beads, α ∈ [1, N_c]. The position of the sth bead in the ith chain is given by r_i(s). The Hamiltonian governing the system behavior is composed by intra- and inter-molecular terms, $Hb$ and $Hnb$⁠, respectively. The intra-molecular interactions that act along the chain’s backbone are given by $βHb=32∑i=1Nc∑s=1Nα−1[ri(s+1)−ri(s)]2/b2$⁠, where b² is the mean squared bond length of an ideal chain, β⁻¹ = k_BT and k_B is the Boltzmann constant. Note that we have assumed identical mean squared bond length and Gaussian configurational statistics for both polymers, but generalizations to more complex situations, including blends of stiff polymers and/or different bond lengths, are straightforward. The inter-molecular interactions are expressed as a functional of the local densities ϕ_A,B(r) according to^23,24

where ρ_o is the bead number density. The first term in the functional represents the repulsion between unlike monomers, quantified by the Flory-Huggins parameter χ. The second term prescribes a finite compressibility for the melt and corresponds to Helfand’s quadratic approximation.²⁵ Parameter κ is related to the inverse isothermal compressibility.^24,26 The field-based description in Eq. (1) can be transformed into a particle-based representation by computing a local density in which a density cloud w(r) is associated with each interacting site of the polymer chains.^23,24,27 This density cloud w(r) decays to zero over a typical distance r_int, which should be large enough such that individual beads interact with a sufficiently large number of particles n_int. In the simulations reported in this work, we have chosen n_int = 14, and we have adopted spherical step-like functions. Specifically, the function w(r) = Cω(r) is defined by ω(r) = 1 if r ≤ r_int/2 and ω(r) = 0 otherwise, r = |r|, and the normalization constant C is defined by the condition that the average density be equal to unity when the chains are ideal. Through these functions, inter-molecular interactions are then expressed as pair-wise interactions, $βHnb=∑i<jU(rij)$⁠, where the sum runs over all polymer segments in the system.^23,24 With this formalism, the intra- and inter-molecular forces acting on each particle can be obtained directly, f^b and f^nb, respectively. As alluded to earlier, the resulting inter-molecular interactions are soft, i.e., the interaction energy remains finite in the limit r → 0, and the local non-crossability constraint is not satisfied and entangled dynamics are not captured.

The topological effects of non-crossability are represented by effective fluctuating interactions, mediated by SS. A schematic representation is given in Figure 1. The SS evolve in pairs, consisting of two rings that move along different chain contours that are connected by a spring. The positions of the rings are denoted by {R_k} and the SS-mediated interactions are included at the same level of description as that of the polymer chains; thus, we use a harmonic interaction between the rings: $βuss=32[Rk−Rnk]2/bss2$⁠, where n_k is the partner of ring k and b_ss characterizes the strength of the spring. The strength of the spring is chosen to be b_ss = b. Our new approach relies on the single-chain model introduced and developed by Schieber et al.^17–19 This model is mathematically well-defined and the associated parameters have a transparent physical meaning. Thus, the SS are only allowed to hop between polymer beads and only one SS occupies a given polymer bead (the so-called “fermionic” case, note that the approach considered in our previous publications does not have this restriction, and more than one SS can occupy a bead, i.e., it is a “bosonic” model). The average number of entanglements per chain, $Z$⁠, is controlled by a pseudo-chemical potential μ which, in turn, is inferred from the molecular weight between entanglements for the particular polymer melt of interest and a specific contour discretization N. Explicitly, using the notation of Schieber,¹⁷ β_e ≡ e^−βμ is chosen to be equal to the average number of beads between entanglements, N_e, i.e., e^−βμ = N_e. To initialize a SS pair, a polymer bead is randomly chosen and the first ring of the SS pair is attached to that segment. The partner ring is then attached to a nearby chain in such a way as to generate a Boltzmann distribution of links.

The conformation of polymers and slip-springs are updated by a hybrid scheme: slip-spring configurations are updated via Monte Carlo (MC) moves every τ_ss time units,³⁵ while polymer configurations are updated via dynamic simulations (Brownian, Dissipative Particle Dynamics (DPD), or any other approach) during this time interval. For the equations governing the temporal evolution of the polymer configurations ${rj}j=1Nbeads$⁠, we have chosen to implement a simple Brownian dynamics approach. Thus, the following stochastic differential equations are integrated:

for j = 1, …, N_beads. In Eq. (2), the forces acting on the beads due to the slip-springs are denoted by f^ss; γ is the corresponding bead friction coefficient. The stochastic variables W_j(t) represent the stochastic influence of the coarse-grained microscopic degrees of freedom and satisfy the fluctuation-dissipation relation. During the MC stage, one performs N_ss(t) + 2N_c attempts to update all SS present at that time t, N_ss(t), as well as to create new SS at chain ends (N_c is number of chains). For each attempt, an existing SS is selected with probability P_SS, or a chain-end is chosen randomly with probability P_ch = 1 − P_SS. Then, one of the following MC moves is attempted:

1. A hopping move involves the transfer of one end of a slip-spring (a ring) from one polymer bead to a neighboring one, along the same polymer contour. We propose to move in one or the other direction along the polymer with equal probability. This move is accepted with probability p = min(1, e^−βΔu_ss).

2. A creation move creates a new slip-spring pair by adding a new ring to a randomly selected unoccupied chain end. The location of the ring’s partner is selected from neighboring polymer beads, according to a Boltzmann distribution, e^−βu_ss. This move is accepted with probability $p=min(1,PSSPchNc2(Nss+1)βe−1)$⁠.

3. A destruction move involves deletion of an existing slip-spring pair if one of the rings reaches the chain end and moves out of the contour with probability 1/2. This move is accepted with probability $p=min(1,PchPSS2NssNcβe)$⁠,

where β_e = e^−βμ. In our implementation, we have chosen $PSS=Z2+Z$⁠. Here, $Z$ is the average number of slip-springs per chain for a given value of the pseudo-chemical potential, and it is obtained from the analytical prediction by Schieber¹⁷ and Khaliullin,¹⁹ namely, $Z=NK+βeβe+1−1$⁠. In the last expression, N_K = N + 1, and N is the total number of polymer beads in a chain. By construction, β_e is chosen to be equal to the average number of beads between entanglements, N_e.

Under this approach, the chain dimensions naturally fluctuate, due the Gaussian chain representation used here, and the SSs favor the movement along the chain’s contour. Two relaxation mechanisms are therefore naturally included in the model, namely, reptation and contour length fluctuations (CLFs). The other main relaxation mechanism, constraint release, arises from the fact that dynamics makes it possible for one ring to leave its chain and, when this occurs, the ring and its corresponding partner can be destroyed. This process mimics the disentanglement at the chain-ends and the constraint release.

It is instructive to highlight the differences between the current implementation and the one used in our previous publications. First, in our previous work, the total number of SS remained constant; the number of slip-springs per chain was allowed to fluctuate, but the ensemble was canonical. Second, SS’s were allowed to move continuously along the polymer backbone, and their evolution was dictated by one-dimensional stochastic differential equations, with an associated friction coefficient, γ_ss. This parameter defined a time scale for the free diffusion of SS’s along the polymer, and it is related to the parameter τ_ss introduced above in our new implementation. An additional difference has to do with the fact that in our previous implementation, a bead could be occupied by more than one SS; in contrast, in the proposed new method, only an individual SS can occupy a bead at a given point in time.

To compare the predictions of the new approach with experimental data in the linear regime, we compute the stress relaxation function, G(t), defined by G(t) = 〈Σ_αβ(t)Σ_αβ(0)〉, where α and β are two orthogonal axes. The stress, Σ_αβ(t), contains contributions arising from intramolecular interactions and slip-spring contributions,

where b_iα(s, t) is the α-component of the bond vector b_i(s) = r_i(s + 1) − r_i(s) that connects the s and s + 1 beads on the ith chain. The time correlation functions required to calculate G(t) are computed using a version of the order-n algorithm.^28,29 From the stress relaxation function, we extract the loss and storage modulus, G′ and G″, by fitting G(t) to a sum of Maxwell modes:³³ $G(t)=G0∑papexp[−t/τp]$⁠, where G₀ = k_BTN_cN/V. Once the parameters a_p and τ_p are obtained, we use the following expressions:

and

It is important to highlight two important physical parameters. The first one is $N¯=(ρoRe3/N)2$⁠, the interdigitation number, which represents the average number of molecules that a chain interacts with and is proportional to the molecular weight of the polymer. This quantity remains invariant independently of what contour discretization N is chosen, and its value is computed from experimental data reported for the polymer of interest. The second parameter, $Re2=(N−1)b2$⁠, denotes the mean squared end-to-end distance of an ideal chain. This parameter sets the unit length scale in our simulations.

We start by performing simulations to corroborate that the implemented grand canonical TIEPOS scheme works correctly and we expect to recover DSM distributions as a self-consistency check. The first aim of our calculations is to compare the simulated $Z$ vs β_e to that predicted by single-chain theories. Simulation results and the corresponding analytical predictions are shown in the main panel of Figure 2. As can be seen, there is good agreement between both. Other characteristics of the system are also described by the theory.^17–19 One can compute the probability distribution of the number of polymer segments between entanglements on a given chain, P(n), the probability distribution P(Q) of the magnitude of connectivity vectors, $Q→$⁠, between SS’s along a given polymer chain, and the probability distribution of the number of entanglements in a given chain, P(Z). The corresponding analytical expressions are^17–19 

where H(x) is the Heaviside step function and δ_i,j is the Kronecker delta. Finally,

where the conditional probability is given by

The results presented in Figures 2 and 3 confirm that our simulations and theory are in good agreement with each other (see supplementary material³⁶ for results obtained using other values of parameter β_e). The consistency between our single-chain and many-chain results serves to underscore the fact that, for homopolymers, we now have the means to jump between different levels of description without sacrificing accuracy.

The predictions of the proposed formalism are next discussed in the context of experimental data for polyisoprene melts by Auhl et al.³⁰ for both the linear and non-linear regimes. The samples considered by these authors were nearly monodisperse and comprised a wide spectrum of molecular weights that encompassed the range from unentangled to well-entangled polymeric liquids. We begin by choosing a reference sample to map experimental parameters to those of our model. We have chosen as reference a sample with molecular weight equal to 33.6 kg/mol (PI-30k). To obtain the corresponding value of $N¯$⁠, we relied on the physical properties reported by Fetters et al. ³¹ for cis-PI at T = 25 °C: from the density (0.91 g/cm³), polymerization index (≃453 monomers), and chain dimension (≃15 nm), we get $N¯=56.18$⁠. We set the end-to-end distance of this polymer as our length scale, R_e = 1, and the contour discretization is chosen to be N = 32. With this choice, the values of ρ_o and b to be used in simulations are specified (see Eq. (1) and intramolecular interaction expression). Different molecular weights can then be represented by choosing proper values of N, with the rest of the parameters fixed. The other important input parameter taken from experimental data is the molecular weight between entanglements, M_e, which is model dependent and determines the height of the plateau modulus. Auhl et al. reported the value of M_e = 4.82 kg/mol, obtained by fitting the tube theory developed by Likhtman and McLeish,¹⁴ to their linear rheology data. However, we found that better agreement between our model and experiment is obtained when the value reported by Fetters et al.³¹ is used, namely, M_e = 3.89 kg/mol.

To study the equilibrium dynamics, the stochastic differential equations (Eq. (2)) are solved using the Euler algorithm. The time is rescaled by the characteristic diffusion time of a system with only random forces, $τ0=l02/D0$ (τ₀ is the diffusion time of a single bead and l₀ is a reference length scale chosen as l₀ = R_e(N = 32) = 1), where D₀ = k_BT/γ is the diffusion constant under these conditions; the rescaled time is given by t^∗ = t/τ₀. In these units, the time step used to integrate the equations is Δt^∗ = 10⁻⁴; the energy is measured in k_BT units. For polyisoprene, we used $N¯=56.18$⁠, χ = 0, and κ = 1.5625. For completeness, all parameters used in our calculations are listed in Table I; the same set of parameters was used for all different molecular weight considered here, both at equilibrium and far from equilibrium.

The only parameter that needs to be fitted is τ_ss, and the value of this time scale is obtained by comparing G(t) of this new implementation and the corresponding function obtained with our previous approach, where the SS friction coefficient was chosen to be γ_ss = γ. The value of τ_ss that gives identical relaxation functions is chosen here. Note that, for completeness, we first performed simulations using our previous implementation for different molecular weights, namely, N = 13, 16, 22, 32, 44, 90. Those with N = 13, 32, and 90 correspond to 13.5, 33.6, and 94.9 kg/mol PI samples reported by Auhl et al., respectively. For PI samples above, the entanglement molecular weight, it was reported that the viscosity scales with molecular weight with exponent 3.6. Our first aim was to corroborate that this scaling is also reproduced by our new model. To this end, the shear viscosity, η, was computed from the stress relaxation function through a Green-Kubo relation,⁶ 

The results of these simulations are shown in the main panel of Figure 4. As can be seen, the 3.6 scaling behavior is recovered by our approach. We have also computed the longest relaxation time τ, defined as $τ=Re2/D$⁠, where D is the diffusion coefficient for the center-of-mass of chains, obtained by the long-time behavior: $limt→∞g3(t)=6Dt$⁠, with g₃(t) = 〈[r_cm(t) − r_cm(0)]²〉, the mean-square center-of-mass displacements. This quantity also displays a scaling behavior with exponent equal to 3.65, consistent with the rheological data (inset in Figure 4).

The loss and storage moduli were determined as indicated in Sec. II. The comparison to experimental data is achieved by scaling the frequency ω and the modulus G₀; the two fitting parameters in the model are determined from the experimental G′ and G″ curves at the point where they intersect each other $(ωc,Gc∗)$⁠. We compute the scaling factors $α=ωcexp/ωcsim$ and $ν=Gc∗,exp/Gc∗,sim$ for the reference sample PI-30k. These scale factors are then applied to all numerical data for the different molecular weights studied here. In Figure 5, simulation predictions and experimental data are presented for the three different samples considered here: PI-14k, PI-30k, and PI-90k (see Table II for details). The average number of entanglements per chain, $Z$ was 3.5, 8.6, and 24.4, respectively. As can be seen, the results from simulations are in quantitative agreement with experimental measurements over a range of time scales that encompasses almost seven decades. As expected, discrepancies start to appear at high frequencies, which correspond to fast processes that cannot be captured by a coarse grained model. Also note that for the PI-90k sample, simulation predictions, although in good agreement with experiments, show slight deviations. We attribute this discrepancy to statistical uncertainty; in order to better resolve the G(t) curve (computed by the time correlation functions of shear stress) at long times, simulations that cover hundreds of relaxation times would be necessary. The noise in G(t) at long times is reflected in the Maxwell modes fit and, therefore, in the complex modulus G^∗(ω).

At this point, we have demonstrated that by fitting only two parameters on a single reference sample, the TIEPOS approach can be used to predict, at a quantitative level, the linear rheological response of other samples with different molecular weights. We now explore the ability of TIEPOS to also describe the non-linear behavior of the samples. To do so, a flow field is added on the differential equations (Eq. (2)) that reproduces the linear velocity profile in shear. The flow proceeds along the x-direction and Lees-Edwards boundary conditions are applied along the gradient direction (y-direction). We compute the transient shear viscosity defined by $η+(t)=σˆxy(t)/γ̇$ for different shear rates $γ̇$⁠. In this expression, $σˆxy(t)$⁠, the shear stress at time t, includes intramolecular, inter-molecular, and slip-spring contributions. Comparison with experiments reported by Auhl et al.³⁰ requires that a time-temperature shift be implemented, G^∗(ω, T) = b_TG^∗(a_Tω, T_o). Their data for η⁺(t) are shifted to a reference temperature of −35 °C. We perform such a transformation by using the horizontal and vertical shift factors, a_T and b_T, obtained from Eqs. (4) and (6) in their report.³⁰ 

Our non-equilibrium simulations focus on the PI-30k sample. We start by computing the linear response from G(t), i.e., $ηlinear+(t)=∫0tG(t′)dt′$⁠. By using the fitted parameters α and ν, as well as the proper shift parameters a_T and b_T, we obtain $ηlinear+(t)$ shown in Figure 6. As expected, this curve is in very good agreement with experimental linear oscillatory shear data. Auhl et al. reported Weissenberg numbers (⁠$Wi=τDγ̇$⁠) associated with their transient shear viscosity measurements. In defining Wi, they used $τD=ηo×Je0$ as the terminal relaxation time, where η_o is the zero-shear viscosity and $Je0$ is the steady-state recoverable compliance, whose values were obtained from measured complex modulus curves (see Ref. 30 for details). As discussed by these authors, the contribution of the different relaxation mechanisms (reptation, CLF, and constraint release) replaces the sharp spectrum of relaxation times with a well-defined dominant reptation time by a continuous spectrum, therefore making the definition of a terminal time difficult. For our shear flow simulations, we have adopted the pragmatic point of view that the terminal time is that associated to the frequency where the G′ and G″ curves intersect each other, $τD,sim−1=ωcsim$⁠. With this value and the reported Wi numbers, we obtain the corresponding simulation shear rates. Again, to make a proper comparison between our predictions and experiments, we need to use the fitted parameters obtained from the linear rheology, α and ν, and the corresponding shift parameters, a_T and b_T. We note here that, as the shear stress as a function of time exhibits pronounced high-frequency fluctuations, we have used an ensemble of 10-20 independent simulations to compute the average curve shear stress vs. time and reduce noise intensity. As can be seen in Figure 6, simulation results are in good agreement with experimental data for all Wi numbers considered here. We emphasize again that no additional parameter adjustments were necessary. The agreement is not perfect, but it is better than that of detailed tube models¹⁴ (compare Figure 6 in this work with Figure 12 on Ref. 30), particularly at high Wi numbers.

The discrepancies that we observe between simulation and experiments could be due to several reasons. First, note that we have used $ωcsim$ to define the terminal time, the latter, as mentioned above, being hard to define. The agreement with data seems to indicate, that the terminal time used is quite close to the “characteristic” terminal time of the real PI sample. Small changes to it will alter the shear rates used in simulations and, therefore, this will be reflected on the computed transient shear viscosity. In general, the simulation predictions are closer to the experimental curves, but not at times before the overshoot. Indeed, by using this grand-canonical approach, the steady-state viscosity is pretty close to that observed experimentally compared with the predictions obtained by using our previous approach (see supplementary material³⁶). The origin of the discrepancy at early times is currently unknown, and studies to clarify this are underway. In Figure 7, we have plotted the average number of slip-springs per chain as a function of time, Z(t), for the different Wi numbers associated to the simulation data in Figure 6. As observed in the single-chain model,^40,41 our simulations also predict a decrease in the number of entanglements as shear rate increases. Interestingly, the time of the maximum of transient viscosity seems to correspond to the time where the rate of slip-spring destruction start to be the largest. At intermediate times, before reaching the steady-state value, Z(t) can be approximated by a power law function, Z(t) ∼ ϑ t^−ε, where ϑ ≈ 10 for all shear rates, and $ε=ε(γ̇)$⁠. The latter has the approximated values {0.25, 0.18, 0.12, 0.075}, from high to low Wi.

We end this section by mentioning that other slip-spring implementations have been proposed in the literature to represent entangled polymers. Chappa et al. originally proposed a similar grand-canonical implementation, as the one adopted in this work, with a fluctuating number of slip-springs, but implemented it on a polymer model consisting of standard DPD beads and finitely extensible nonlinear elastic (FENE) springs both for slip-springs and polymer bonds. As such, the correspondence to a more coarse-grained representation (and the corresponding prescription for the pseudo-chemical potential for the slip-springs) was not established.^9,34 Müller-Plathe et al.¹³ and Panagiotopoulos et al.³² implemented methods with a constant number of slip-springs, where the temporal evolution of the material relied on alternating Dissipative Particle Dynamics and slip-spring Monte Carlo move sequences (each having a fixed number of steps). Such implementations used a relatively large number of MC steps to update the slip-springs distribution, thereby allowing the SS configuration to relax to a “deeper” minimum while polymers remained fixed in space. Such implementations could lead to a collapse of the molecules and might not be applicable far from equilibrium. A quantitative comparison with experimental data has not been presented for all such approaches, and it is difficult to anticipate at this point how each model will perform in the context of realistic polymeric materials.

We have introduced a new TIEPOS approach for quantitative descriptions of the rheological behavior of polymer melts, both in the linear and non-linear regimes. In contrast to a previous formalism introduced by our group, in this work we rely on a grand-canonical representation that allows the number of slip-springs (entanglements) to fluctuate. The results of the new model are shown to be in agreement with a single-chain theory, thereby facilitating the selection of some of the parameters that are necessary to describe a given system. The results are also in agreement with those of our previous, canonical implementation, except at high shear rates, where the number of slip springs decreases. We have shown that by fitting only two scaling factors to data for a single sample, it is possible to predict the rheology of samples having different molecular weights, both at equilibrium and under shear, without any further parameter adjustments. The results presented here were focused on a set of well-characterized linear polyisoprene melts. In previous work, we showed that TIEPOS provides a good, quantitative description of polystyrene samples; in that case, experimental data were taken from different laboratories, and therefore each sample under consideration was fit. Both, monodisperse and bidisperse systems were considered and agreement with experiment was highly encouraging. All these results suggest that the main relaxation mechanisms that arise in entangled polymers are correctly captured by the proposed TIEPOS approach. Furthermore, the fact that TIEPOS and DSM give nearly identical results is encouraging and suggests that both approaches can be integrated to perform systematic coarse-graining. DSM calculations are simple to code and fast. The TIEPOS formalism is more demanding, but because it relies on a many-molecule representation of the system, it can go beyond single-chain theory and be used to describe inhomogeneous materials, such as composites or multiblock polymers. A systematic study of TIEPOS for such systems will be presented in the near future. Moving forward, it will also be important to develop and test approaches that allow one to go coarse-grain atomistic representations into TIEPOS model and DSM.

This work was supported by U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Science and Engineering Division. An award of computer time was provided by the INCITE program of the Argonne Leadership Computing Facility. We gratefully acknowledge the computing resources provided on Blues, a high-performance computing cluster operated by the Laboratory Computing Resource Center at Argonne National Laboratory. J.D.S. acknowledges the support of the National Science Foundation, NSF-CBET Fluid Dynamics No. 1438700, and the Army Research Office, No. W911NF-11-2-0018. Marat Andreev gratefully acknowledges the support from NIST through a CHiMaD postdoctoral award.