Primitive-path statistics of entangled polymers: mapping multi-chain simulations onto single-chain mean-field models

We present a method to map the full equilibrium distribution of the primitive-path (PP) length, obtained from multi-chain simulations of polymer melts, onto a single-chain mean-field ‘target’ model. Most previous works used the Doi–Edwards tube model as a target. However, the average number of monomers per PP segment, obtained from multi-chain PP networks, has consistently shown a discrepancy of a factor of two with respect to tube-model estimates. Part of the problem is that the tube model neglects fluctuations in the lengths of PP segments, the number of entanglements per chain and the distribution of monomers among PP segments, while all these fluctuations are observed in multi-chain simulations. Here we use a recently proposed slip-link model, which includes fluctuations in all these variables as well as in the spatial positions of the entanglements. This turns out to be essential to obtain qualitative and quantitative agreement with the equilibrium PP-length distribution obtained from multi-chain simulations. By fitting this distribution, we are able to determine two of the three parameters of the model, which govern its equilibrium properties. This mapping is executed for four different linear polymers and for different molecular weights. The two parameters are found to depend on chemistry, but not on molecular weight. The model predicts a constant plateau modulus minus a correction inversely proportional to molecular weight. The value for well-entangled chains, with the parameters determined ab initio, lies in the range of experimental data for the materials investigated.

includes fluctuations in all these variables as well as in the spatial positions of the entanglements. This turns out to be essential to obtain qualitative and quantitative agreement with the equilibrium PP-length distribution obtained from multi-chain simulations. By fitting this distribution, we are able to determine two of the three parameters of the model, which govern its equilibrium properties. This mapping is executed for four different linear polymers and for different molecular weights. The two parameters are found to depend on chemistry, but not on molecular weight. The model predicts a constant plateau modulus minus a correction inversely proportional to molecular weight. The value for wellentangled chains, with the parameters determined ab initio, lies in the range of experimental data for the materials investigated.

Introduction
Many dynamic properties of concentrated polymer solutions and polymer melts are determined by the topological constraints that the macromolecular chains impose on one another, called entanglements. Since each entanglement involves at least two chains, it is not obvious a priori whether these properties can be predicted from a single-chain mean-field level of description. In early attempts, Edwards [1] and de Gennes [2] introduced reptation theory, which was adopted by Doi and Edwards [3][4][5][6] as the basis of their nonlinear flow model. In reptation theory, the discrete entanglements of a probe chain with other chains are replaced by a continuous 'tube'. In the classical models of de Gennes and Doi and Edwards, equilibrium dynamics are governed by one-dimensional diffusion (reptation) of the probe chain along the tube axis. Reptation may cause a chain end to move inward and vacate its tube segment 9 , upon which that segment is destroyed, or to move outward and explore new constraints, creating a new, randomly oriented tube segment. The extended reptation model for linear viscoelasticity of linear homopolymers, developed by Likhtman and McLeish [7], includes the additional relaxation mechanisms of contourlength fluctuations, constraint release and monomer density equilibration. Further mechanisms have been postulated to obtain agreement with linear viscoelastic data for more complex molecular architectures [8,9] or polydisperse systems [10,11]. Also many modifications of the Doi-Edwards model have appeared to improve nonlinear flow predictions [12][13][14][15][16][17][18][19][20]. These models involve various assumptions and approximations, and are mostly able to predict data for one type of experiment (i.e. a specific kind of material and/or a specific flow). We agree with Likhtman [21] that this situation is undesirable and that efforts should instead be directed toward a multi-scale approach, in which reptation models are derived, with clear assumptions, from more-detailed single-chain or multi-chain models. In a recent paper [22], such a coarsegraining procedure was applied to the mobile slip-link model (MSM) proposed by Schieber and Horio [23]. This resulted in a generalization of the chain free energy for the reptation model of Read et al [20]. Another example is the derivation of a continuum model for unentangled polymers from a multi-chain finitely extensible nonlinear elastic (FENE) dumbbell model [24]. For a review of multi-scale methods for polymers, see [25].
In the present work, again in the spirit of the multi-scale approach [21], we show that detailed multi-chain simulations can be mapped onto the MSM to determine the parameters that govern the equilibrium statistics of this model. This leaves only one adjustable parameter, a characteristic time scale, in the MSM for flow predictions. The mapping is done on the level of the primitive-path (PP) network. Algorithms have been developed to reduce each chain in a multi-chain system simultaneously to a set of straight line segments, which connect the two chain ends through successive entanglements [26][27][28][29][30][31]. These straight line segments define the PP of the chain, which should correspond to the tube axis in reptation models or a set of entanglement connector vectors in slip-link models. The original PP algorithm of Everaers et al [26] minimizes an energy function under the constraints of fixed chain ends and mutual uncrossability of chains. Since the number of entanglements and their distribution along the chain are not produced as output, the average number of monomers between entanglements can only be estimated with the help of a model. Later algorithms, such as Z and Z1 by Kröger and co-workers [28,30] and CReTA by Tzoumanekas and Theodorou [31], minimize the contour length of each chain (subject to the above constraints) through geometric operations. These geometric algorithms keep track of the monomers that take part in entanglements, and therefore the average number of monomers between entanglements can be obtained in a modelindependent way.
Several groups have attempted to obtain equilibrium properties of the Doi-Edwards tube model and its successors from multi-chain PP networks. Most of them focused on the average molecular weight between entanglements M e [26,[31][32][33][34][35][36]. This quantity is not predicted by the tube model and cannot be measured experimentally. However, by inserting the Doi-Edwards stress tensor in the Green-Kubo relation, one obtains the short-time modulus as G(0) = ρ RT /M e . Here ρ is the density, R is the universal gas constant and T is the absolute temperature. It is now widely accepted among reptation theorists that the plateau modulus G 0 N , measured in experiments, is lower than the short-time Green-Kubo result by a factor 4 5 . 10 The average molecular weight between entanglements can thus be estimated from the experimental plateau modulus as M e = 4 5 ρ RT /G 0 N . The corresponding average number of monomers between entanglements is where M m is the molar mass of a monomer. The common procedure in tube-based PP analysis is to compare N G e to the result of another calculation, which takes as input the average squared end-to-end distance R 2 ee eq and the average PP length L pp eq at equilibrium, obtained from the PP network: (2) 10 Likhtman and McLeish [7] derived this prefactor based on the assumption that the monomer density rapidly equilibrates to minimize tension differences between tube segments.
Here N m = M/M m is the number of monomers per chain. Both N G e and N pp e depend on the assumptions that the tube axis (the PP) is a random walk, that each chain has the same number of tube segments Z , and that each tube segment has the same number of monomers N e (or number of Kuhn steps N E ). The calculation of N pp e additionally depends on the assumption that each tube segment has the same length, the 'tube diameter' a = √ N E a K , where a K is the Kuhn length. The effects of some of these assumptions on the result obtained for N pp e are discussed in appendix A.
By contrast, in the single-chain mean-field slip-link model introduced by Schieber [37], entanglements are described as discrete objects (slip-links) whose number fluctuates, and is governed by a chemical potential bath in a grand-canonical ensemble. Also the numbers of Kuhn steps in the strands as well as the orientations and lengths of the entanglement connector vectors fluctuate. The equilibrium distributions of all these stochastic variables can be obtained analytically [37]. Foteinopoulou et al [32] showed that they agree very well with the distributions in multi-chain united-atom polyethylene systems, which were equilibrated by a connectivity-altering Monte Carlo algorithm and subsequently reduced to PP networks. In similar work on polyethylene and cis-1,4-polybutadiene, Tzoumanekas and Theodorou [31] found that the distribution p eq (N ) of the number of Kuhn steps per strand deviates from the sliplink model prediction only at small numbers, N / N eq 0.5. They were able to account for this difference by introducing an excluded-volume interaction between entanglements on the same chain, which is not present in the slip-link model. The resulting distribution was verified for poly(ethylene terephthalate) [33] and polystyrene [34] melts; hence it appears to be universal.
In an early version of the slip-link model, called the continuous slip-link model [38][39][40], the numbers of Kuhn steps in the strands could take noninteger values. Sliding dynamics (SD) of Kuhn steps through entanglements was driven by chemical potential differences between strands instead of tension differences, as is generally postulated in reptation models 11 . In the more rigorous version, called the discrete slip-link model [41], the numbers of Kuhn steps in the strands are treated as integers. This allows a mathematically rigorous formulation of creation and destruction of entanglements at the chain ends by SD, in terms of discrete jump processes [41]. Moreover, creation and destruction of entanglements everywhere along the chain (constraint dynamics) is incorporated in a self-consistent way [41]. With these two relaxation mechanisms, based on minimal assumptions, the slip-link model has been shown to predict linear viscoelasticity and dielectric relaxation of monodisperse linear, bidisperse linear and branched polymer melts [41][42][43] as well as linear viscoelasticity of crosslinked polymer networks [44] quantitatively 12 .
Khaliullin and Schieber [45] derived analytic expressions for the first two moments of the equilibrium PP-length distribution of the continuous slip-link model. The chain-length dependence of the first moment, L pp eq , was found to agree qualitatively with results obtained by Foteinopoulou et al [32] from PP analysis of united-atom polyethylene systems. Quantitative agreement with these multi-chain results was obtained for L pp eq as a function of chain length in the discrete slip-link model, which Khaliullin and Schieber calculated numerically. However, for the continuous as well as the discrete slip-link model, the second moment L 2 pp eq turned out 11 The formulation in terms of chemical potentials naturally prevents complete disentanglement, without introducing a hypothetical force at the chain ends as proposed by Doi and Edwards [4]. 12 The only exception is a discrepancy with dielectric relaxation data for highly dilute long chains in a matrix of short chains, which has not been explained by any entanglement theory yet [43].
to be larger than in the multi-chain systems. Nevertheless, it is encouraging that the calculations were done with a value of the entanglement-density parameter β which was consistent with the multi-chain results, and that this yielded the prediction G(0) = 3.0 MPa for the short-time modulus of polyethylene at T = 413 K, in reasonable agreement with the experimental plateau modulus G 0 N = 2.6 MPa [46,47]. Note that the short-time modulus G(t = 0) of a model may be larger than the plateau modulus G 0 N if some rapid relaxation processes are present in the level of description of that model. We believe that the success of the slip-link model is due to the fact that it accounts for many fluctuations, which are known to exist in multi-chain systems [31,32]. From this point of view, it is surprising that the PP analysis based on the tube model, in which the only quantity with a finite-width distribution is the segmental orientation, has shown reasonable consistency between the monomer density estimates N pp e and N G e for several chemistries [26,[31][32][33][34][35][36]. However, N pp e consistently disagrees with the average number of monomers per PP segment, determined directly from PP-network data (from the number of kinks, presumably corresponding to entanglements, per PP) without using any model assumptions [28, 31-36, 48, 49]. This average is commonly estimated according tō as proposed by Kröger [28]. The expressionN ≈ (0.4-0.5)N pp e ≈ (0.4 − 0.5)N G e captures the results from many groups [50].
Hoy et al [48] introduced new definitions of N pp e andN , which they called the 'M-estimators' M-coil and M-kink, respectively. The improvement over the original 'S estimators' S-coil, given by (2), and S-kink, given by (3), which were used by other groups [26,28,[31][32][33][34][35][36]49], is in the rate at which they converge to a plateau when plotted against chain length. In [48, figure 5] it is shown that the plateau values of the M estimators are the same as those of the S estimators. These new definitions do not solve, but helped to highlight, the discrepancy between N pp e andN , which is probably due to the neglect of fluctuations in the tube model underlying the S-coil definition and the M-coil definition of N pp e . Support for this explanation is given in appendix A.
Thus the tube-based PP analysis exhibits self-consistency between model-dependent estimates of the monomer density, using data either from multi-chain PP networks or from plateau modulus measurements, but inconsistency between these estimates and modelindependent results derived from multi-chain PP networks. In view of this problem, we regard the slip-link-based PP analysis by Khaliullin and Schieber [45], which is not only selfconsistent, but also to a large extent in agreement with multi-chain simulations, as a more solid foundation to build on. One problem is that the slip-link model predicts a PP whose length distribution is broader than that obtained from multi-chain simulations. Another problem is that this PP prediction has random-walk statistics, whereas it has been known for some time that correlations between PP segments exist in multi-chain systems [28,31]. In other words, the PP has a persistence length larger than the average entanglement spacing. It is clear that some physics, currently unaccounted for, needs to be incorporated to solve these problems.
Entanglement spatial fluctuations (ESFs) are a candidate for the missing physics in singlechain entangled polymer melt models. Spatial fluctuations of network nodes have a long history in crosslinked network models [51][52][53][54][55]. The same idea has been applied to entangled melts, for example in the reptation model of Read et al [20], who showed that the mean PP behaves like a random walk with an added harmonic bending modulus. Qin et al [56] extended this model with a semiflexibility correction to contour-length fluctuations, which improved agreement with molecular dynamics simulations in terms of the tangent correlation function along the PP. ESFs were also included, on more detailed levels of description, in the slip-spring simulation of Likhtman [57] and the MSM of Schieber and Horio [23]. Here we use a recently proposed version of the MSM, in which the ESFs become anisotropic in flow [58]. We refer to the original model of Schieber and Horio [23] as the isotropic MSM and to the new model as the anisotropic MSM. At equilibrium, where the statistics of the two models are equivalent, we simply use the term MSM. However, it is important to note that the dynamics of ESFs affects the form of the stress tensor and hence the relaxation modulus. It is shown elsewhere that a specific time evolution of ESFs is required to satisfy four consistency criteria, which arise naturally in singlechain mean-field models [58]. We come back to this issue in section 4. The slip-link model without ESFs [37][38][39][40][41][42][43][44]59] is called the fixed slip-link model (FSM). Here we use equilibrium statistics, based on the thermodynamics of the MSM, to analyze PP network data from multichain simulations.

Mobile slip-link model
The MSM has been formulated on two different levels of description, both of which are singlechain and mean-field [22,23]. On the more-detailed level, the set of stochastic variables is where Z is the number of strands of the chain, N i is the number of Kuhn steps in the ith strand, Q i is the end-to-end vector of the ith strand and X i is the end-to-end vector of the ith virtual spring. Each virtual spring symbolizes the confinement potential acting on an entanglement; X i points from the center of the potential (the anchor) to the entanglement, as shown in figure 1. Schieber et al [58] added a tensorial variable n i for each entanglement, which describes the size and the directionality of the ESFs, i.e. n −1 i describes the strength and shape of the confinement potential. At equilibrium, all ESFs are isotropic (n i = nδ ∀i) and therefore these tensors are left out of the set of variables here. The equilibrium probability density is then given by the modified Maxwell-Boltzmann expression where J is a normalization constant, µ e is the chemical potential conjugate to the number of entanglements Z − 1 and is the Gaussian chain free energy. The virtual number of Kuhn steps n is one of the three adjustable parameters of the MSM. The second is the activity Together, n and β determine all static properties of the model. One more adjustable parameter, the time scale τ K , is needed for the dynamics [38][39][40][41][42][43]59]. The total number of Kuhn steps N K is not adjustable, since we know how to calculate it for a given chain chemistry; see for example appendix B. The total number of Kuhn steps is preserved by the Kronecker delta in (5). Compared to the isotropic MSM proposed by Schieber and Horio [23], we have incorporated not only entanglement anisotropy [58] but also, following the work of Andreev et al [59], dangling-end contributions to the free energy. Like the entangled strands, the dangling ends are described as random walks. This means that their average length is zero, but higher moments of their end-to-end vector distribution are nonzero. The dangling ends then contribute to the stress, and consequently to the relaxation modulus. As shown in section 3.1, these contributions have no significant influence on the plateau modulus predicted for wellentangled chains. Also shear flow predictions are unaffected [59]. However, keeping track of the conformations { Q 1 , Q Z } of the dangling ends might be important in extensional flow at high strain rates [59].
If the PP is defined as the set of instantaneous strand end-to-end vectors { Q i }, then its length has the equilibrium distribution We use the discrete model unless noted otherwise. For the continuous model, the sums over Kuhn-step numbers are replaced by integrals and the sum over Z runs from 1 to ∞. Since (8) is identical to the result obtained for the FSM, the discrepancy with multi-chain simulations [45] remains. However, in the MSM, the definition of the PP is not unique. Other possible definitions are the 'anchor path', formed by the set of anchor connector vectors {q i }, or the 'mean path', formed by the set of strand end-to-end vectors {Q i } that minimizes the free energy for a given number of strands, Kuhn-step numbers and anchor connectors [22]. To derive the equilibrium length distributions of these alternative definitions of the PP, we integrate p eq ( ) over all possible conformations of the virtual springs {X i }, while keeping the variablesω fixed [22,23]. This is then the set of variables of the MSM on the coarse-grained level of description, where the actual positions of the entanglements are no longer known. We can change variables tô by using the relation [22] where for Z > 1 and for Z = 1. The two variable setsω andˆ thus correspond to the same level of description. The integration over all possible virtual-spring conformations leads to the equilibrium probability density with the coarse-grained chain free energy where | . . . | represents taking the determinant. The equilibrium distribution of the anchor-path length is then The analogous quantities p eq (ˆ ), F(ˆ ) and pQ eq (L pp ) are obtained by making the replacementŝ in (14)- (16). It is important to realize that both the anchor path and the mean path have a persistence length greater than the average strand length. The anchor path has the longest persistence length, since all elements of the matrix σ −1 in (15) are nonzero, which means that all anchor connector vectors are directly coupled. The matrixσ −1 is tridiagonal, which means that each mean strand end-to-end vector is directly coupled to its nearest neighbors only. Therefore we anticipate that the statistics of the mean path are closer to those of the PP in multi-chain systems, since correlations between PP segments have been observed in these systems [28], which decay exponentially with increasing separation between the segments [31]. The short-time modulus of the anisotropic MSM, in both of the coarse-grained variable spacesω andˆ , is [58] Here n c is the number density of chains. For n = 0, we recover the result G FSM (0) = n c k B T Z eq obtained for the discrete FSM [41] with dangling-end contributions to the stress. The discrete isotropic MSM [23] yields the same result, hence it has a front factor f = G iso (0)/G FSM (0) = 1. According to (12) and (19),σ decreases with increasing n, i.e. with increasing size of ESFs. This means that the anisotropic MSM has a front factor in accordance with the multi-chain entangled melt simulations of Masubuchi et al [60] and with crosslinked network models [51][52][53][54][55].

Results
The derivation of analytic expressions for the distributions p q eq (L pp ) and pQ eq (L pp ) or for their moments is complicated by a combination of two factors. The first is the coupling of the strand conformations by the full matrix σ −1 in (15) or by the tridiagonal matrixσ −1 in F(ˆ ), the chain free energy expressed in the mean-path variables, which is obtained from (15) by the replacements (17)- (19). The second complicating factor is the coupling of the Kuhn-step numbers by the Kronecker delta in (14) or in its analogue p eq (ˆ ), the equilibrium probability Table 1. Properties of monodisperse cis-1,4-polyisoprene multi-chain systems. Multiple independent simulation ensembles were used for each system and N c is the combined number of chains in these ensembles. Calculations with this code show that (for n > 0) the distribution of the anchor-path length, p q eq , is broader than that of the instantaneous PP length, p Q eq , which in turn is broader than the distribution of the mean-path length, pQ eq . This suggests that the mean path might be closer to the PP obtained from multi-chain simulations than either the anchor path or the path defined by the instantaneous strand end-to-end vectors, since Khaliullin and Schieber [45] showed that the equilibrium PP-length distribution from multi-chain simulations is narrower than p Q eq . In the next few subsections, we show that this is indeed the case. Quantitative agreement is obtained between the cumulative distribution (22) and its counterpart P mc eq (L pp ) from PP analysis of multi-chain systems, with a unique set of model parameters {n, β}. These parameters depend on the chain chemistry, but not on the molecular weight.

Monodisperse cis-1,4-polyisoprene
Li et al [61] performed molecular dynamics simulations of monodisperse cis-1,4-polyisoprene systems with different chain lengths. They used a coarse-grained description in which each monomer (containing four backbone carbon atoms) is represented by one superatom. The PP network was obtained by the geometric algorithm Z1 [30]. Table 1 lists the properties of the systems, which are labelled as in [61]: the number after 'PI' is the number of beads per chain, which is equal to N m + 1. The number of Kuhn steps per chain N K is calculated using the number  (2), i.e. the S-coil definition [26,48], using equations (4)-(6) of Li et al [61] with parameter values from their table 2. Circles:N according to (3), i.e. the S-kink definition [28,48], using data from as derived analytically for the discrete (either fixed or mobile) slip-link model [23,37], is included in table 1 for all systems. The value of β is determined by fitting PQ eq (L pp ) to P mc eq (L pp ), as explained below.
It appears to have gone unnoticed that the results of Li et al [61] seem to disagree with those of other groups. The solid line in figure 2 shows N pp e versus N m according to the S-coil definition [48], given by (2). We calculated this as the ratio of two functions, which Li et al used to fit the 'tube diameter' a pp (N m ) = R 2 ee eq / L pp eq and the average PP length per monomer b pp (N m ) = L pp eq /N m for systems with different N m . These functions are given in [61, equations (5) and (6)]. The circles in figure 2 showN according to the S-kink definition [48], given by (3), which we calculated using data extracted from [61, figure 9(d)]. The ratio between the two estimates is roughly two, in accordance with previous work [50]. They must converge to the same plateau values as the M-coil and M-kink estimators [48], which were reported as N pp e = 185 andN = 76 by Li et al [61]. The interesting result, however, is that the plateau value of N pp e is greater than N G e (the dashed line in figure 2, obtained from (1) using the plateau modulus measured by Abdel-Goad et al [62]) 13   at least 50% 14 . If this result were correct, then this would suggest that the agreement between N pp e and N G e in previous works [26,[31][32][33][34][35][36]50] might be a coincidence. This would imply that the tube model is not fully self-consistent, since N pp e and N G e are both model-dependent. Curiously, Li et al [61] foundN according to the M-kink definition, which is modelindependent, to be close to N G e . Figure 3 presents the cumulative mean-path distribution PQ eq (L pp ) from discrete MSM calculations, using an ensemble size of 4 × 10 4 chains, and the cumulative PP-length distribution P mc eq (L pp ) from multi-chain simulations, for systems PI40-PI400. The same is shown for PI500-PI1000 in figure 4. For PI260-PI1000, multi-chain simulations were performed in addition to those originally published. The systems were equilibrated by the method used in [61], which consists of running a molecular dynamics simulation, combined with a connectivityaltering Monte Carlo algorithm, and monitoring the decay of the autocorrelation function of the end-to-end unit vector. The Monte-Carlo algorithm accelerates the equilibration process and also increases the ensemble size by generating independent configurations along the way. Our new results correspond to longer equilibration times and larger ensemble sizes; see [61, table 1] and our table 1. Convergence of the PP-length distribution was confirmed for PI800 and PI1000. The MSM results shown in figures 3 and 4 were obtained with n = 2.4 and β = 15.7, which were determined by manual fitting of the PP-length distributions of systems PI260-PI1000. Most of the multi-chain results up to PI1000 are described accurately by the MSM with this parameter set. Significant discrepancies appear for PI40 and PI80, where Z eq 3. For these systems, 14 As shown in the supplementary data (available from stacks.iop.org/NJP/16/015027/mmedia), for the simulations of Li et al [61], the M-coil estimator does not seem to be converged even for the system with the longest chains. This means that the plateau value of N pp e might be even larger than reported in their paper.  N K , which is about 10% smaller than the FSM prediction, lim N K /β→∞ L pp eq = √ 2/(3β) = 0.206 from (A.10). The variance is about a factor four lower than the FSM prediction, lim N K /β→∞ L 2 pp eq − L pp 2 eq = N K a 2 K from (A.11). Figures 5 and 6 include the average and the variance of p mc eq (L pp ) from the original results of Li et al [61]. For N K 200, these results deviate from the linear dependence on N K , predicted by the MSM. This highlights the critical importance of equilibration in PP analysis, which was already discussed in detail by Hoy and Robbins [63]. We expect that improved equilibration of systems PI1500 and PI2000 will similarly bring their PP-length distributions closer to the MSM predictions. Note also that the new multi-chain simulations improve the molecular-weight dependence of N pp e from the S-coil definition, as shown by the triangles in figure 2. Since the average Kuhn-step number is β + 1 for well-entangled chains [37] and N m /N K = 2.21 as derived in appendix B, the value β = 15.7, obtained by fitting the MSM to multi-chain PP-network data, corresponds to a monomer density N e = 37.0. This is nearly half of the tubemodel estimate N G e = 72.7; hence it is close to what one sees in earlier work [31][32][33][34][35][36]50]. The short-time modulus, calculated using (20), is plotted versus 1/N K in figure 7, in the dimensionless form We find that the expression describes the MSM data accurately for N K β, i.e. for Z eq 2. The inset of figure 7 zooms in on the well-entangled regime (N K β) to show the close agreement between (24) and (25). The same behavior is observed for atactic polystyrene and cis-1,4-polybutadiene, which are discussed in sections 3.2 and 3.3, respectively. Note that (25) does not reduce to (24) in the limit n → 0. Therefore we expect the agreement to break down at some lower value of n. We have verified, both theoretically and numerically, that the exact result (24) has the correct asymptotic behavior. Figure 7 also includes the analytic result obtained for the discrete FSM [41] with dangling-end contributions to the stress. The front factor f of the MSM is given by (21). We divide (25) by (26) and take the limit N K /β → ∞, which leads to as an estimate of the asymptotic front factor for well-entangled chains. With the parameters n = 2.4 and β = 15.7, we find f ∞ = 0.54 for cis-1,4-polyisoprene. This is consistent with results from multi-chain entangled melt simulations [60] and close to the theoretical front factor for crosslinked networks with a functionality of four: f = 1/2 [51][52][53][54][55].
Around the same chain length where (25) deviates from the exact result (24), the contribution of unentangled chains becomes important. This is evidenced by the fact that G(0) approaches the analytic result for a system made up entirely of such chains, which are equivalent to Gaussian dumbbells, It is interesting to look at the short-time modulus of the anisotropic MSM with the danglingend conformations coarse-grained out, similar to previous work on the FSM [37][38][39][40][41][42][43] and the isotropic MSM [23]. In the limit N K /β → ∞, the fraction of dangling ends goes to zero and therefore G(0) is the same as when the chain ends are tracked. As N K decreases, the fraction of dangling ends increases; G(0) then drops because dangling ends do not contribute to the stress. This is shown in figure 7 for cis-1,4-polyisoprene. Similar results have been obtained for the other materials. Since the number of Kuhn steps in a strand (in the discrete slip-link model) is at least one, G(0) is exactly zero for N K 2. Liu et al [64] observed a molecularweight dependence in the plateau modulus derived from neutron spin-echo (NSE) measurements on polybutadiene [65], which looks qualitatively similar to our results without dangling-end contributions. However, the derivation of Liu et al makes use of a fit of the tube model to the static structure factor obtained from NSE. Since the tube model neglects many of the fluctuations present in NSE, whereas the MSM captures more (but maybe not all) of them, no definitive conclusions can be drawn at this point. The dimensionless short-time modulus for well-entangled chains, lim N K /β→∞ M K G(0)/ (ρ RT ) = 0.0325, lies within the range of experimental plateau modulus data reported by Fetters et al: M K G 0 N /(ρ RT ) = 0.0387 in [46], M K G 0 N /(ρ RT ) = 0.0320 in [47] and M K G 0 N /(ρ RT ) = 0.0244 in [62]. In the last reference, the plateau modulus is essentially constant down to a (weight-averaged) molecular weight of 13.5 kDa, which corresponds to 1/N K = 0.011. This is still in the regime where we observe only a weak chain-length dependence of G(0). Therefore the MSM agrees with these experimental data.

Monodisperse atactic polystyrene
We present multi-chain PP-network data for monodisperse atactic polystyrene here, polydisperse cis-1,4-polybutadiene in section 3.3, and finally monodisperse and polydisperse polyethylene in section 3.4. These systems were all equilibrated by a connectivity-altering Monte Carlo algorithm and reduced to PP networks by means of the CReTA algorithm [31]. The simulations were performed at the temperatures and with the ensemble sizes given in table 2. The packing length was estimated for all these systems, and was found to be in very good agreement with published data [66].
In the polystyrene simulation, each CH 2 -CH(C 6  chain is N b = 2N + 2 = 4002. For more details concerning the atomistic structure and coarsegraining of atactic polystyrene, we refer to Spyriouni et al [34]. With the characteristic ratio C N = 9.8, extracted from the equilibrated chain conformations, we derive the Kuhn-step length a K = 18.1 Å and the number of Kuhn steps per chain N K = 282. Figure 8 shows that the mean-path distribution of the discrete MSM is practically identical to the PP-length distribution obtained from the multi-chain simulations, consistent with our results for cis-1,4-polyisoprene in section 3.1. The parameter n for polystyrene is smaller than that found for cis-1,4-polyisoprene: n = 1.0 versus n = 2.4, which means that polystyrene exhibits weaker ESFs. One might speculate that this is related to the difference in chain stiffness (C ∞ = 5.3 for cis-1,4-polyisoprene [67]) and/or the presence of bulky side groups such as the aromatic rings in polystyrene. Applying our PP analysis method to a wider range of materials could provide more insight into the influence of chemical properties on n. The ratio n/(β + 1) is of the same order for the two materials: n/(β + 1) = 0.098 for polystyrene and n/(β + 1) = 0.14 for cis-1,4-polyisoprene. (In the limit N K /β → ∞, β + 1 is the average number of Kuhn steps per strand [37,41].) With the model parameters n = 1.0 and β = 9.2 from our PP analysis, (24) yields the dimensionless short-time modulus M K G(0)/(ρ RT ) = 0.0628 for atactic polystyrene with N K = 275. Using M K = 727 Da [59] and p-V-T data from the literature [68], we obtain G(0) = 327 kPa at T = 175 • C. This seems to agree with the dynamic modulus measurements of Schweizer et al [69] on polystyrene with a molecular weight M = 200 kDa, which corresponds to N K = 275, at the same temperature. However, the experimental plateau modulus is hard to locate exactly; a more stringent test is to run a linear viscoelastic simulation with the MSM, using the model parameters determined ab initio, and compare the results to experimental data. Andreev et al [59] fitted the dynamic modulus data of Schweizer et al [69] with the discrete FSM including dangling-end contributions to stress, from which they obtained β = 15.1. Inserting this value in the analytic expression (26) for the dimensionless short-time modulus of their  (27), we obtain the front factor f ∞ = 0.61 for well-entangled chains, which is slightly higher than that found for cis-1,4-polyisoprene.

Polydisperse cis-1,4-polybutadiene
We now consider a polydisperse cis-1,4-polybutadiene system. In the simulation, each CH x group was coarse-grained to a united atom. The distribution of chain lengths is uniform with a minimum of 600 and a maximum of 1400 united atoms. The corresponding number-averaged molecular weight and polydispersity index areM n = 13.5 kDa andM w /M n = 1.05, respectively. More details can be found in the papers of Tzoumanekas and Theodorou [31] and Gestoso et al [70]. The characteristic ratio C 1000 = 4.6 is almost identical to C ∞ = 4.7, as reported in [70]. Using C 1000 , we find the Kuhn-step length a K = 8.8 Å and the number of Kuhn steps N K ∈ [78,182]. The calculation is analogous to that presented in appendix B for cis-1,4polyisoprene.
In a modified version of our code, we select the N K of each chain according to the known molecular weight distribution of the multi-chain system, and subsequently generate the chain conformations in the same manner as before, so that the ensemble has the equilibrium statistics of the discrete MSM. In section 3.1, the parameters n and β were found to be independent of molecular weight. This means that we should be able to fit the PP-length distribution, obtained from the polydisperse multi-chain simulation, with a single parameter set. Figure 9 shows that this is indeed the case. Consistent with our results for monodisperse melts in sections 3.1 and 3.2, the mean-path distribution of the discrete MSM is nearly indistinguishable from the multi-chain PP-length distribution. The parameters found are n = 2.0 and β = 12.0. Their ratio n/(β + 1) = 0.15 is very close to that of cis-1,4-polyisoprene.   [46,47]. The front factor is identical to that of cis-1,4-polyisoprene: f ∞ = 0.54, again in agreement with results from multi-chain entangled melt simulations [60] and crosslinked network models [51][52][53][54][55].

Monodisperse and polydisperse polyethylene
In the polyethylene simulations, each CH 2 group was coarse-grained to a united atom. We first consider a monodisperse system, for which the ensemble size is 32 000 chains as shown in table 2. The number of united atoms per chain is N = 1000 and the number of bonds per chain is N b = N − 1 = 999. Details can be found in the work of Tzoumanekas and Theodorou [31]. From the characteristic ratio C ∞ = 7.8, which is slightly lower than the C ∞ = 8.27 reported previously [32], we calculate the Kuhn-step length a K = 14.6 Å and the number of Kuhn steps per chain N K = 87. Figure 10 shows excellent agreement between the PP-length distributions from the multichain simulations and our fit using the MSM. However, we obtain n/(β + 1) = 1.6. A ratio greater than one means that the characteristic time scale of the virtual-spring conformations is longer than that of the strand conformations. Therefore our coarse-graining procedure, in which we integrate out the virtual-spring conformations from the detailed MSM [22,23], might not be appropriate for polyethylene. The low value of β, which implies a low average number of Kuhn steps per strand, also raises some issues. Most importantly, the assumption of random-walk statistics for the strands becomes unrealistic. We would probably need to replace the Gaussian chain free energy F( ), given in (6), by a semiflexible free energy. To avoid confusion, we stress here that the semiflexibility of the PP [28,31] is a separate issue. This property of the PP arises naturally in the MSM when coarse-graining the Gaussian chain free energy F( ) to the mean-path free energy F(ˆ ), as explained in section 2. Surprisingly, the parameters obtained from our PP analysis of polyethylene suggest that the chain, not just its PP, is semiflexible.
Semiflexible polymer models generally treat the PP as a continuous space curve, but also assume a uniform Kuhn-step (or monomer) density. Our results and those of Khaliullin and Schieber [45] show that fluctuations in the Kuhn-step numbers are important for PP statistics. Thus the development of a semiflexible model with nonuniform monomer density would be useful for PP analysis of polyethylene, and possibly other polymers as well. Another problem for β ∼ O(1) is that it is no longer safe to approximate the Kuhn-step numbers as integers. This approximation is naturally avoided when using a semiflexible model. As seen in figure 10, merely switching to the continuous MSM does not solve our problem. The result is very close to that of the discrete MSM, but now n/β = 3.8. (For the continuous slip-link model, the average number of Kuhn steps per strand in the limit N K /β → ∞ is β [37].) Using the model parameters {n, β} = {3.2, 1.0} from our PP analysis to calculate the short-time modulus, and extrapolating to infinite degree of entanglement as explained above, we find lim N K /β→∞ G(0) = 5.68 MPa for polyethylene at T = 413 K. This is higher than the experimental value of G 0 N = 2.6 MPa, obtained at the same temperature by Fetters et al [47], although they did not report the molecular weight. With N K = 87, as in the multi-chain simulations considered here, (20)  are 384 000 chains for the system with number-averaged molecular weightM n = 14 kDa and 40 000 chains for the system withM n = 84 kDa. Simulation details for the latter system can be found in the work of Uhlherr et al [71]. Figures 11 and 12 also include the MSM predictions, using the same parameters as for the monodisperse polyethylene melt, which are listed in table 2. The agreement is reasonable, which again implies that these parameters are independent of molecular weight.

Discussion
PP analysis has received considerable attention [26, 31-36, 48, 49] as a tool to map multi-chain simulations onto single-chain mean-field models. The first step in PP analysis is to choose a target: a single-chain model to map onto. The common choice is the Doi-Edwards tube model. However, this model neglects several fluctuations, whose importance in multi-chain systems has been known since the early days of PP analysis [31,32]. Later work based on the tube model has suffered from a discrepancy between model-dependent estimates of the monomer density and direct estimates based solely on data from multi-chain PP networks. 'Fluctuations' are sometimes pinpointed as the source of this discrepancy, but explanations are typically based on models that do not include all the fluctuations observed in multi-chain simulations [50]. As shown in appendix A, including more of these fluctuations leads to a decrease of model-based estimates of the monomer density, which brings them closer to direct estimates.
Khaliullin and Schieber [45] used the FSM for PP analysis. The FSM contains fluctuations in the orientations of entanglement strands (which are the only fluctuations present in the tube model) and additionally in the lengths of the strands, the number of entanglements per chain and the number of monomers or Kuhn steps in each strand. With a single adjustable parameter, β, the FSM describes the average PP length at equilibrium as a function of molecular weight, as extracted from multi-chain simulations, and predicts a short-time modulus consistent with experiments. However, the shape of the PP-length distribution disagrees with multi-chain results [45].
In this paper, we have presented a PP-analysis method based on the MSM [22,23,58]. In addition to the fluctuations already present in the FSM, the MSM includes ESFs. The strength of these fluctuations is governed by one extra adjustable parameter, n. The MSM allows different possible definitions of the PP. Our results show that the path taken by the strand end-to-end vectors when they are in their mean positions, which minimize the free energy, corresponds to the PP obtained from multi-chain simulations. The shape of the equilibrium PP-length distribution can be described quantitatively by adjusting n and β. These parameters do not depend on the molecular weight, but do depend on the chain chemistry; see table 2. The value of β obtained is smaller than the one found by Khaliullin and Schieber [45] for the FSM. This means that ESFs decrease the short-time modulus, in agreement with results from multi-chain simulations [60]. We investigated four linear polymers: cis-1,4-polyisoprene, atactic polystyrene, cis-1,4-polybutadiene and polyethylene. The results for polyisoprene, polystyrene and polybutadiene might suggest a universal ratio n/(β + 1) ≈ 0.1, while for polyethylene n > β and β is of order one, which indicates that a semiflexible chain model would be more appropriate for this material. The fact that the parameters n and β are similar for the structurally similar polyisoprene and polybutadiene (see table 2) is encouraging. The short-time moduli, predicted by the MSM with the parameter values from our PP analysis, agree with plateau modulus measurements for well-entangled chains. Finally, in the MSM, the mean strand end-to-end vectors are correlated: each strand is directly coupled with its nearest neighbors. This is at least qualitatively consistent with the correlations between PP segments, observed in multi-chain systems [28], which decay exponentially along the PP [31]. The same holds for the reptation model of Read et al [20], but that model does not contain fluctuations in the degree of entanglement or the monomer density. It is shown in the supplementary data (available from stacks.iop.org/NJP/16/015027/mmedia) that the PP-length distribution predicted by the model of Read et al does not agree with multi-chain simulation results. Also in the supplementary data, the joint distribution of PP length and end-to-end length is examined for the Doi-Edwards model, the FSM, the MSM and a multi-chain FENE simulation. The results demonstrate that the MSM is very close to the multi-chain simulation.
We conclude that the MSM is an excellent target model for PP analysis. To the best of our knowledge, the literature contains no precedent of a self-consistent mapping of multichain simulation data onto a single-chain mean-field model with the level of qualitative and quantitative agreement demonstrated here. The determination of the model parameters n and β can be done at a small computational cost by the code included in the supplementary data. Moreover, reduction of multi-chain systems to PP networks is nowadays a standard procedure. A direction for future work is the determination of the one remaining adjustable parameter in the MSM, the characteristic time τ K for sliding of Kuhn steps through entanglements. Hence ab initio rheology predictions may now almost be within reach.
The main message from this paper is that all fluctuations, contained in the MSM, are necessary to obtain agreement with the PP statistics from multi-chain simulations. A more comprehensive investigation of the importance of ESFs is presented in a separate paper involving two of us [58], where the following four criteria for single-chain mean-field entangled polymer melt models are considered: (i) Consistency with a multi-chain level of description. Masubuchi et al [60] compared multichain simulations in which the spatial positions of entanglements were allowed either to fluctuate or to be displaced affinely with the macroscopic flow. Their results show that the plateau modulus decreases when ESFs are included. Based on studies of crosslinked network models [51][52][53][54][55], it may be expected that the correct implementation of ESFs in single-chain models would result in a front factor f < 1 in the relation between the shorttime modulus and the average molecular weight between entanglements. This would lead to a lower average molecular weight between entanglements predicted by these models, in accordance with the results of PP analysis. (ii) Consistency with nonequilibrium thermodynamics. Several formalisms have been proposed for the development of models that do not violate the laws of thermodynamics. The most successful nonequilibrium thermodynamics formalism is GENERIC [72][73][74][75][76], which, among other things, restricts the form of the stress tensor. (iii) Consistency with the stress-optical rule. It is well known from experiments on entangled polymers that, under circumstances where the chain conformations are Gaussian, the stress tensor is linearly proportional to the refractive index tensor [77,78]. This is known as the stress-optical rule. (iv) Consistency between the Green-Kubo relation and the response to small perturbations. The relaxation modulus G(t) can be calculated from the autocorrelation function of the stress tensor using the Green-Kubo relation. Linear response theory requires the autocorrelation of the stress tensor at equilibrium to be consistent with the stress increase due to an infinitesimal step strain imposed on the system.
The FSM satisfies the last three of these criteria. It can be cast in the GENERIC structure and the resulting stress tensor satisfies the stress-optical rule [79]. The short-time modulus G(0) from the Green-Kubo relation is given in [41]. The proof that this agrees with the modulus after an infinitesimal step strain is straightforward. However, the discrete slip-link model does not include ESFs and thus violates criterion (i).
Schieber and Horio [23] showed that the incorporation of isotropic ESFs in the slip-link model does not affect the short-time modulus. Even worse, the stress tensor has a contribution from the ESFs, which violates the stress-optical rule. Thus the isotropic MSM of Schieber and Horio satisfies only two of the four criteria. However, by casting the MSM in the GENERIC form, it can be shown that the stress tensor contribution from the ESFs disappears if these fluctuations become anisotropic in flow according to a specific time evolution [80]. This time evolution is a generalization of an idea proposed for crosslinked polymers by Ronca and Allegra [81]. With the generalized Ronca-Allegra dynamics, the short-time modulus decreases when the size of ESFs increases [58]. The Green-Kubo relation remains consistent with the response to small perturbations; hence all four criteria are satisfied [58]. The relation of our PP analysis to other recently published work is discussed in appendix C.
We start with a single-chain model in which the number of strands Z is constant and each strand has the same constant number of Kuhn steps N E . The orientations and, contrary to tube models, also the lengths of the strands are allowed to fluctuate. The equilibrium probability density for the chain conformation is then After integrating (A.1) over the two angles in a spherical coordinate system, the distribution of From the former, we derive the average strand length and, with Z = N K /N E , the average PP length L pp eq = Z Q eq = 8 3π Substitution of R 2 ee eq = N K a 2 K leads to Multiplying both sides of this expression by the number of monomers per Kuhn step, which is a constant determined by the chain chemistry, we find This is smaller than the result for the tube model by a factor 8/(3π) ≈ 0.85, which is due to the fact that we do not neglect strand-length fluctuations. The fact that N e moves closer to the monomer densityN , obtained directly from multi-chain simulations (see the introduction), gives us hope that incorporating more of the fluctuations that are present in these simulations [32] will lead to further improvement. If Kuhn steps are allowed to slide through entanglements, then entanglements can be destroyed (if the last Kuhn step in a dangling end slides through its entanglement) or created (if a Kuhn step at the very end of a chain gets entangled with another chain). Therefore any selfconsistent model that includes fluctuations in the number of entanglements should also include fluctuations in the numbers of Kuhn steps in the strands. The slip-link model incorporates both of these in addition to fluctuations in the orientations and lengths of the strands, which were already present in the model described above. The equilibrium probability density for the FSM is where J is a normalization constant and is the Helmholtz free energy of the chain [37].
For the continuous FSM, in which the Kuhn-step numbers can take noninteger values, the equilibrium PP-length distribution is The mth moment of this distribution is given by In order to solve this equation analytically, it is convenient to start by integrating over the strand lengths, then use the procedure from [23, appendix 2] to integrate over the Kuhn-step numbers, and finally sum over the number of strands. The first two moments are [45] L pp eq and Hence the parameter β, which for the continuous slip-link model corresponds to the average number of Kuhn steps per entanglement, β = N K /( Z eq − 1), is given by  .4), we see that the numerical prefactor decreases further when (in addition to fluctuations in the strand end-to-end vectors) fluctuations in the entanglement number and the Kuhn-step numbers are accounted for.

Appendix B. Kuhn-chain parameters of cis-1,4-polyisoprene
For sufficiently long chains, the average squared end-to-end distance at equilibrium is where C ∞ is the characteristic ratio, N b is the number of bonds in the chain backbone and l 2 is the mean squared bond length [82]. The Kuhn chain is defined as a random walk with the same contour length as the real chain, i.e.
where N i m and L i m are the number of monomers of type i and the contour length per monomer of type i, respectively.
The chemical structure of 1,4-polyisoprene is With the molar mass of a monomer M m = 68.119 g mol −1 , we calculate the molar mass of a Kuhn step M K = 151 g mol −1 . The ratio R 2 ee eq /M = a 2 K /M K = 0.677 Å 2 mol g −1 agrees with the results of Li et al [61] from simulations, R 2 ee eq /M = 0.665 Å 2 mol g −1 , and Fetters et al [46] from experiments, R 2 ee eq /M = 0.679 Å 2 mol g −1 .

Appendix C. Relation to recent work
• Uchida et al [86] showed that PP-network data for many flexible and semiflexible polymer melts and solutions can be characterized by scaling relations based on the packing length [46,47,66], suggesting that the PP length is a static property.
• Likhtman et al [87] obtained the PP-length distribution for a toy model consisting of a single Rouse chain on a cubic lattice, where the grid lines act as rigid obstacles. From this they calculated N pp e and compared it to estimates from several procedures based on dynamics, in which the monomer density appears as a parameter: fitting the relaxation modulus G(t) with an expression derived from slip-spring simulations [57,88], fitting the mean squared displacement of the middle monomer with tube-model scaling laws in different regimes, and fitting the dynamic structure factor in NSE measurements with an empirical formula [65]. The estimated monomer densities were consistent neither with each other, nor with N pp e .
• The works of Uchida et al [86] and Likhtman et al [87] are compatible with our MSMbased PP analysis, which shows that the average number of Kuhn steps per strand (i.e. the parameter β) can be obtained from a purely static property of the multi-chain system: the equilibrium PP-length distribution. Determination of the characteristic time scale τ K for sliding of a Kuhn step through an entanglement, which is the single adjustable parameter for the dynamics of the MSM [38][39][40][41][42][43][44]59], is beyond the scope of the present paper. A few works have been devoted to the determination of dynamic properties of reptation models, such as the tube segment survival probability [89][90][91]. We argue that the MSM is a more suitable target to map dynamic multi-chain properties onto than the tube model, again because the MSM includes more of the fluctuations present in multi-chain systems.
• In a detailed analysis of united-atom polyethylene simulations, Anogiannakis et al [92] showed that a large fraction of all interactions between chains have very short lifetimes. Looking at a snapshot of a PP network, many of these weak links might be counted as entanglements, although their effects on the PP or the stress are probably negligible. Anogiannakis et al proposed a method to select the most persistent links, whose lifetimes exceed some critical value. These should correspond to the entanglements in single-chain mean-field models. In our approach, we do not need to impose any such selection criterion.
The value of β obtained from fitting the PP-length distribution corresponds only to those interchain interactions that affect the PP statistics. In following papers, we will verify whether the same β also leads to accurate predictions of linear viscoelasticity and nonlinear flow.
• Andreev et al [59] investigated the three approximations made in the discrete FSM, which might be inappropriate for nonlinear flow predictions: neglect of convective constraint release, instantaneous relaxation of dangling ends upon destruction by SD, and unbounded extensibility of the strands. They incorporated convective constraint release by means of a chain-pairing algorithm [93,94], finite deformation and relaxation of dangling ends through Brownian dynamics of the chain ends, and finite extensibility by using Cohen's [95] Padé approximation to the inverse Langevin function for the strand free energy. It was found that shear flow predictions are insensitive to all three approximations, and agree quantitatively with experimental data. The same is true for uniaxial elongation up to a Hencky strain between two and three. After that, the discrete FSM predicts a viscosity overshoot and severe thinning, which are not observed experimentally.
• Andreev et al [59] concluded that some physics are missing in the discrete FSM-and by extension, in all less-detailed single-chain mean-field entangled polymer models. We now know that the anisotropic MSM satisfies four consistency criteria for such models [58], of which one was violated by the original discrete FSM (see the discussion at the end of section 4). The present paper shows that the anisotropic MSM predicts the right PP statistics at equilibrium. It remains to be seen whether it also fixes the discrepancy with nonlinear extensional flow predictions.