# Static and dynamic properties of large polymer melts in equilibrium

###### Abstract

We present a detailed study of the static and dynamic behavior of long semiflexible polymer chains in a melt. Starting from previously obtained fully equilibrated high molecular weight polymer melts [Zhang et al. ACS Macro Lett. 3, 198 (2014)] we investigate their static and dynamic scaling behavior as predicted by theory. We find that for semiflexible chains in a melt, results of the mean square internal distance, the probability distributions of the end-to-end distance, and the chain structure factor are well described by theoretical predictions for ideal chains. We examine the motion of monomers and chains by molecular dynamics simulations using the ESPResSo++ package. The scaling predictions of the mean squared displacement of inner monomers, center of mass, and relations between them based on the Rouse and the reptation theory are verified, and related characteristic relaxation times are determined. Finally we give evidence that the entanglement length as determined by a primitive path analysis (PPA) predicts a plateau modulus, , consistent with stresses obtained from the Green-Kubo relation. These comprehensively characterized equilibrium structures, which offer a good compromise between flexibility, small , computational efficiency, and small deviations from ideality provide ideal starting states for future non-equilibrium studies.

## I Introduction

A fundamental property of polymer melts containing long linear chains is that they are entangled. As the stiffness of chains increases, the entanglement effect becomes stronger, i.e. the entanglement length is shorter. Complex topological constraints in polymer melts play an essential role for dynamic, and rheological properties. For studying such properties and phenomena in an out-of-equilibrium state it is important to begin with a well characterized equilibrium ‘sample’ of very long polymer chains in a melt. It is the purpose of this study to provide this.

According to Flory’s argument, the excluded volume interactions become screened de Gennes (1979); Yamakawa (1971) when the concentration of polymer solutions exceed the chain overlap concentration. Therefore polymer chains in a melt eventually behave statistically as ideal chains, as if excluded volume effect would no longer be important. However, Wittmer and his co-workers Wittmer et al. (2007, 2011) have pointed out that there are noticeable deviations from an ideal chain behavior due to the incompressibility constraint of the melt. For fully flexible polymer chains in a melt based on lattice and continuum models, bond fluctuation model (BFM) and bead-spring model, respectively, such deviations are indeed seen. This finding is confirmed by a recent Monte Carlo study of polymer melts using BFM in Ref. Hsu (2014a) while the deviations are less visible as the chain stiffness starts to play a role for polymers. Therefore, we provide a detailed study of the conformational properties of long bead-spring polymer chains in a melt as the chain stiffness is taken into account, where we especially study to what extent polymer chains behave as ideal chains.

It is well known that for short unentangled chains in a melt, the motion of monomers can be approximately described by the Rouse model de Gennes (1979); Doi and Edwards (1986); Kalathi et al. (2014, 2015). If the polymer chains become long enough such that the effects of entanglements start to become important, movements of chains at the intermediate time and length scales are confined to a tube-like region, created by surrounding chains and depending on the corresponding entanglement length . The dynamic behavior within this time frame is well described by the tube model of de Gennes, Doi and Edwards de Gennes (1979); Doi and Edwards (1986); Rubinstein and Colby (2003). Each polymer chain is assumed to move back and forth (reptation) along the contour of an imaginary tube around the so called primitive path. Although ample evidence of reptation scaling predictions is given by previous Monte Carlo and molecular dynamics simulations Baumgärtner and Binder (1981); Kremer and Grest (1990); Paul et al. (1991); Wittmer et al. (1992); k. Kremer and Grest (1992); Kopf et al. (1997), a complete picture still is lacking. This is mostly due to the limitations of available equilibrated systems of huge chain length and the long relaxation times covering several orders of magnitude.

Recently, the authors of Ref. Zhang et al. (2014) developed a novel and very efficient methodology for equilibrating high molecular weight polymer melts through a sequential backmapping of a soft-sphere coarse-grained model Vettorel et al. (2010); Zhang et al. (2013) from low resolution to high resolution, and finally the application of molecular dynamics (MD) simulations of the underlying bead-spring model (see Appendix). Therefore, a further investigation of the static and dynamic scaling behavior predicted by theories de Gennes (1979); Yamakawa (1971); Doi and Edwards (1986) for huge systems in the highly entangled regime has become easily accessible. Therefore, the aim of this paper is to give a deeper understanding of static and dynamic behavior of large semiflexible polymer chains in a melt, and compare our numerical results whenever it is possible to theoretical predictions in the literature. We mainly focus on polymer melt system containing semiflexible polymer chains of sizes , , with the Flory characteristic ratio . The chains are modelled as standard bead-spring chains with a bond bending interactions parameter . For details of the model we refer to the appendix. All results quoted refer to chains with a bending constant of unless otherwise noted. Having such big polymer melt systems at hand we have the possibility to analyze the linear viscoelasticity as characterized by the stress relaxation modulus, and estimate the entanglement length from the standard expression of the plateau modulus . It is also interesting to check whether is equivalent to the estimate of through the primitive path analysis (PPA) Everaers et al. (2004).

The outline of the paper is as follows: Sec. II describes the static conformational structures of polymer chains in a melt, and compares them to those for ideal chains. Sec. III describes the motions of polymer chains in a melt at different characteristic time scales, and verifies the scaling laws predicted by the Rouse model, and the reptation theory de Gennes (1979); Doi and Edwards (1986); Rubinstein and Colby (2003). The detailed structure investigation of the primitive path of chains is given in Sec. IV. Studies of linear viscoelasticity of polymer melts are given in Sec. V. Finally, our conclusions are summarized in Sec. VI.

## Ii Static properties of equilibrated polymer melts

Let us first look at the estimates of the mean square end-to-end distance and the mean square radius of gyration given by

(1) |

and

(2) | |||

where is the position of monomer of chain number while is the center of mass (c.m.) of the -th polymer chain in a melt, and the average includes an averaging over all independent equilibrated configurations. Results of and plotted versus are shown in Fig. 1 for polymer melts containing chains of sizes , , and . Here the root-mean square bond length . We see that , and as one would expect for ideal chains.

The conformational behavior of individual polymer chains of size in a melt can also be described by the probability distributions of end-to-end distance and gyration radius , and , respectively. For ideal chains where , the probability distribution of is a Gaussian distribution,

(3) |

Although there exists an exact theoretical prediction Yamakawa (1971); Fujita and Norisuye (1970); Denton and Schmidt (2002) for the probability distribution of it is much more complicated to evaluate. However, it has been checked Vettorel et al. (2010); Hsu (2014a, b) that the formula suggested by Lhuillier Lhuillier (1988) for polymer chains under good solvent conditions in -dimensions is still a good approximation for ideal chains (), i.e.,

(4) |

where and are (non-universal) constants, and the exponents and are linked to the space dimension and the Flory exponent by and . Here is the des Cloizeaux exponent des Cloizeaux (1975) for the osmotic pressure of a semidilute polymer solution, and is the Fisher exponent Fisher (1966) characterizing the end-to-end distance distribution.

The probability distribution of any observable is normally obtained numerically by accumulating the histogram over all configurations and all chains of size , and then normalizing the histogram such that

(5) |

In Fig. 2, we present the normalized probability distribution () as a function of () for polymer melts of three different chain sizes , , . Note that an angular average over all directions has been included in . We see the nice data collapse for both and , and they are described very well by the following two -independent normalized distribution functions obtained from Eqs. (3), (4), and with shown in Fig. 1,

(6) | |||

and

(7) | |||

where the parameters , , and the normalization factor are determined numerically by a least-squares fit.

For understanding the connectivity and correlation between monomers the conformations of linear chains of contour length in a melt are usually described by the average mean square internal distance, ,

(8) |

where is the chemical distance between the monomer and the monomer along the identical chain. It is generally believed that the theoretical prediction of mean square internal distance for polymer melts consisting of semiflexible chains in the absence of excluded volume effect described by a freely rotating chain model is Flory (1969)

(9) |

with

(10) |

In the limit , the bond-bond orientational correlation function therefore decays exponentially as a function of chemical distance between any two bonds along a linear chain Grosbeg and Khokhlov (1994); Rubinstein and Colby (2003),

(11) |

where is the so-called persistence length which can be extracted from the initial decay of .

As , Eq. (9) gives the asymptotic behavior of the mean square end-to-end distance of a FRC equivalent to the behavior of a freely jointed chain

(12) | |||||

(13) |

where is so-called Flory’s characteristic ratio Flory (1969), and is the Kuhn length.

Results of scaled by (), obtained by taking the average over independent polymer melts containing chains to reduce fluctuations at large , are shown in Fig. 3. For , we see the nice data collapse for chains of different sizes . The universal scaling behavior for is nearly in perfect agreement with the theoretical prediction of for semiflexible chains in the absence of excluded volume effect described by a freely rotating chain (FRC) model. However, a slight deviation from the predicted curve for FRC occurs for . This deviation becomes more prominent as the flexibility of polymer chains increases due to the correlation hole effects that the correlation hole is deeper for more flexible chains. Note that here we do not take the bond-bond orientational correlation between two successive bond vectors, in Eq. (9), as a fitting parameter Auhl et al. (2003), but rather we estimate directly from the equilibrated configurations of polymer melts.

The correlations between two bonds along an identical chain at a chemical distance for , , are shown in Fig. 4a. As it was clarified in Refs. Wittmer et al. (2007); Hsu et al. (2010); Hsu (2014b), the asymptotic decay of as a function of for dense melts and at the point is not a single exponential as predicted by Eq. (11), but rather a power law decay, for , due to excluded volume effects. Therefore only the initial decay of is meaningful for the estimation of the persistence length . However, the crossover point shifts to larger value of as the chain stiffness increases, i.e. the range over which the exponential decay holds extends. We also check how the profiles of probability distribution of bond angles vary with increasing chain stiffness. Using Eq. (5), is estimated by accumulating normalized histograms of between two successive bonds along a chain. We see that in Fig. 4b, the distributions have a bimodal form. For fully flexible chains () in a melt, there exists one peak occurring at due to the competition between the excluded volume effect and the flexibility. As the chain stiffness increases ( increases), a second peak starts to develop at , and the position where the peak is located shifts to a smaller value of . For an ideal chain in a dilute solution, one should expect that

(14) |

This is also shown in Fig. 4b for comparison.

The scattering from single chains in a melt in equilibrium is shown in Fig. 5. In Fig. 5a we see that for small (, ) in the Guinier regime, then a crossover occurs to the power law of Gaussian coils (ideal chains), with for . Here for using Eq. (11). Though our chains are moderately stiff () the short range initial rigid-rod regime for is hardly visible, thus allowing them still to be taken as a model for flexible polymers. In order to clarify whether single chains in a melt behave as ideal chains we show the structure factors in a Kratky-plot in Fig. 5b. The Debye function de Gennes (1979); Cloizeaux and Jannink (1990); Schäfer (1999); Higgins and Benoit (1994) describing the scattering from Gaussian chains,

(15) |

is also presented in Fig. 5b for comparison. The deviations from ideality are clearly recognized near for rather flexible chains (, ) of size , and a minimum value is reached in the Kratky-plot as increases, in agreement with the previous work Wittmer et al. (2007); Hsu (2014b). As a first conclusion one can state that polymer melts of chains with a stiffness parameter offer a good compromise for modeling highly flexible polymers while at the same time minimizing deviations from ideality, which significantly impair the use of simple models for fully flexible chains.

## Iii Dynamic properties of equilibrated polymer melts

The dynamic behavior of polymer chains in a melt or solution is usually characterized by the mean square displacement (MSD) of monomers. The theoretical predictions of the dynamic scaling behavior of MSD given by the reptation theory de Gennes (1979); Doi and Edwards (1986) show that the crossover behavior occurs at different time scales, the characteristic time , the entanglement time , the Rouse time , and the disentanglement time (in the ideal case where the chain length is very large). However, all simulations and experiments support due to the reason that contour length fluctuation, constrains release and correlation hole effects shift the crossover to the asymptotic behaviors to very long chains Doi (1983); Milner and McLeish (1998); McLeish (2002); Likhtman and McLeish (2002).

Three quantities describing the dynamic properties of polymer chains in a melt are listed as follows: the mean square displacement of a monomer,

(16) |

the mean square displacement of monomers with respect to the corresponding center of mass (c.m.),

and the mean square displacement of the center of mass

(18) |

Note that in Eq. (16) only half of the monomers in the middle of each chain are considered in order to suppress the fluctuations caused by chain ends Kremer and Grest (1990); k. Kremer and Grest (1992), while all monomers in each chain are considered in the calculation of the center of mass {Eq. (18)}. The corresponding scaling predictions of , , and are given by Kremer and Grest (1990); Pütz et al. (2000)

Our extensive molecular dynamics results of , , and up to for polymer chains of sizes , in a melt are shown in Fig. 6. The best fits of the theoretical predictions given in Eq. (LABEL:eq-g123-s) are shown by solid lines for comparison. The characteristic time scales where is the LJ time unit (see Appendix), , and for are determined by the intersection points of two lines from the results of in Fig. 6a. They correspond to the crossover points between two scaling regimes are pointed out by arrows also in Fig. 6bcd. The disentanglement time is determined from the intersection between the fitting straight lines of for and for , respectively, since we should expect that for . The characteristic time estimated from for is which is compatible with the direct measurement. If we estimate the entanglement length from characteristic time scales , , and determined by the scaling predictions of the mean square displacement for (Fig. 6a and 6d), we get and if we assume that . Both estimates are consistent with results from PPA and from the relaxation plateau modulus within error bars (see Table 1). The two estimates are deviating by about from the expected value . If we fit our data with , we get which is underestimate. Thus our data perfectly fit experiments that and show the limitations of the asymptotic theory.

According to the theoretical predictions, we see that in Fig. 6a, at . At , (assuming that a Rouse chain of monomers is relaxed) Kremer and Grest (1990); k. Kremer and Grest (1992), where the entanglement effect starts to set in and monomers in an identical chain are restricted to move only along the contour of an imaginary tube of diameter and contour length until reaching for . Since the tube itself is a random walk with a step length , the displacement of a monomer at is thus . In the case of , we find that our data of () follow the power law about three decades for , a much longer time window than observed so far via simulation. For , the polymer chain slides back and forth along the tube-like regime and results in a second regime which is predicted by the reptation theory Doi and Edwards (1986). After reaching the disentanglement time (reptation time) , a chain has moved a distance comparable to its own size for (see Fig. 6d). The initial tube is completely destroyed and another new tube-like regime will appear depending on the polymer chain size or polymer molecular weight. Finally monomers diffuse such that for .

Results of , , and (Fig. 6abc) show that they are all independent of for . Furthermore in that regime. For , either the size is still too short or the statistics for long relaxation time are insufficient, the expected scaling law is only seen slightly, while for , for is barely reached. However, such a proof for or even longer chain lengths might only be possible with further improved soft and hardware Anderson et al. (2008); Glaser et al. (2015).

## Iv Comparison between the original chain conformations and the primitive path

In order to understand the structural differences between the original path and the primitive path (pp) of polymer chains in a melt, we implement the same primitive path analysis proposed by Everaers et al. Everaers et al. (2004) based on the concept of Edwards’ tube model Edwards (1967) to identify the primitive path of each polymer chain in a melt Sukumaran et al. (2005); Kröger (2005); Shanbhag and Larson (2005); Tzoumanekas and Theodorou (2006); Hoy et al. (2009); Everaers (2012). A detailed discussion regarding to self-entanglements, local self-knot effect, and finite-size effect is given in Ref. Sukumaran et al. (2005); Moreira et al. (2015).

Since the motion of a chain is confined in a tube-like regime with fluctuation due to entanglements with other chains (see Sec. III), the primitive path of the chain is the contracted contour of an imaginary tube without any other chain crossing when all endpoints are fixed in space. In this analysis, topologies of chains are kept and chains are assumed to behave as random walks along their primitive paths. The mean square end-to-end distance of chains therefore remains the same as that for the original paths of chains, i.e., , and

(20) |

Here is the Kuhn length, is the contour length, and is the average bond length of the primitive path. The so-called entanglement length defined by the number of monomers per Kuhn segment of the primitive path is then

(21) |

Quantitatively, the primitive paths of all polymer chains in a melt are determined by slowly cooling the system toward and minimizing the energy of the system Everaers et al. (2004); Moreira et al. (2015). In the simulation, two ends of chains are fixed and the intrachain excluded volume interactions as well as the bond bending interactions are switched off while the interchain interactions are kept. In the case where the intrachain excluded volume is kept, Sukumaran et. al. Sukumaran et al. (2005) have found that the difference of the estimate of between these two cases is within error bars. Results of the bond-bond orientational correlation function , and the normalized histogram of bond angles , for the primitive paths of polymer chains in a melt with are shown in Fig. 7. The initial decay of described by an exponential decay up to is shown by a dashed line with . Since the endpoints of chains are fixed, without considering the interchain interactions and thermal fluctuations, chains are stretched out when the bond springs try to reduce the average bond length from to . This effect is stronger at the short length scale () where the result of show some deviations from the fitting curve if we take a closer look. The stretching conformations of chains are also observed from the normalized histogram of bond angles shown in Fig. 7. The distribution of still has a bimodal form, but the range of shrinks from (Fig. 4b) for the original paths to for the primitive paths in the case of . The distance between two peaks decreases as decreases.

Results of the mean square internal distance for the original and the primitive paths of polymer chains in a melt with are shown in Fig. 8. Since the endpoints of chains are fixed, one should expect that results of for both paths approach to the same value with increasing . It is indeed seen in Fig. 8. If we use where with in Eq. (9), we see that results of for the primitive path can still be well described by the FRC. We also check the distributions of bond length {Eq. (20)} for the primitive paths and show that the distribution is simply a normal (Gaussian) distribution of () given by

(22) |

where is the standard deviation of , and is the mean value of (Fig. 9). The distributions of the entanglement length , , for , , and , and for are shown in Fig. 9b. We see that does not depend on . The position of the peak of (fig. 9c) corresponds to the estimate Moreira et al. (2015) of . Results of through the PPA are listed in Table 1 for three different chain sizes and for .

## V Viscoelasticity

The viscoelasticity of polymer melts is normally characterized by the stress relaxation modulus as a function of relaxation time . For , since the dynamics of chains can be described by the Rouse model while reaches a plateau value depending on the entanglement length, or the molecular weight between entanglements predicted by the reptation theory Doi and Edwards (1986); Doi (1980) for where chains are assumed to move in a tube-like regime due to entanglements. Finally, entangled chains are relaxed for and starts to deviate from the plateau.

In order to clarify whether the entanglement length estimated from stresses using the standard expression of the plateau modulus is equivalent to determined through PPA mentioned in Sec. IV, we perform MD simulations to estimate the stress relaxation modulus . Two methods are considered here. One is from the stress autocorrelation function (SAF) of off-diagonal elements of the preaveraged stress tensor for fully equilibrated polymer melts Lee and Kremer (2009); Lee et al. (2010). The components of stress tensor taking the pairwise potential and the three-body potential into account are defined via the virial theorem:

(23) | |||||

where and are the mass and the th component of the velocity vector of the th bead, respectively, and is the th component of the force vector acting on the th bead by the th bead. Using the Green-Kubo relationship Sen et al. (2005), the stress relaxation modulus

(24) |

where the off-diagonal element . In order to reduce the strong noise in SAF Lee and Kremer (2009); Lee et al. (2010), is defined by

(25) |

where the preaveraged stress tensor

(26) |

In our simulations, we choose MD steps with the time step .

(plateau) | ||
---|---|---|

The other method is to measure the normal stress decay after deforming polymer chains in a melt by a small step strain elongation, since linear viscoelastic properties are associated with near equilibrium measurements of the system where the configurations of polymer chains are not moved far away from their equilibrium states. In our simulations, this is done by applying cycles of uniaxial elongation to deform the simulation box with a strain rate (holding each chain in a tube-like regime) at each cycle such that at the end the simulation box is elongated in the -direction (), but shrunk in the -, -directions (). Here the volume of the simulation box is kept fixed, , and the stretch ratio with such that the system is in the linear viscoelastic regime. Using the stress-strain formulas for classical rubber elasticity Treloar (1986), the stress relaxation modulus

(27) |

Results of scaled by with estimated by the plateau value of are shown in Fig. 10. The estimates of are also listed in Table 1. They are in perfect agreement with the estimates through PPA within error bars. In Fig. 10a, is estimated from Eqs. (24)-(26) for polymer melts consisting of chains of sizes , , and , and for . Due to the difference between microscopic structures of independent equilibrated polymer melts, we observe that the plots of as a function of show slightly different scenarios for different sets of data (not shown). Therefore, besides taking the preaverage of for the estimate of , we shall also take the average of over independent sets of data although our systems are quite large. For , the scaling law predicted by the Rouse model is verified. As increases, the curves of for three different sizes first reach a plateau for , then start to deviate from it depending on the chain size . Since , the range over which extends with increasing . However, in Fig. 10b, we only focus on the case of and compare the results of obtained from two different measurements, Eqs (24) and (27). For the second measurement, two values of the strain rate are chosen, and . We see that only depends on for . For , results of estimated from the normal stress tensor are consistent with the estimates from .

## Vi Conclusion

In this paper, we have studied bead-spring chains in a melt at a monomer density by extensive molecular dynamics simulations using the ESPResSo++ package Halverson et al. (2013). We investigate the static and dynamic properties of polymer chains in a melt. For fully equilibrated large polymer melts, we observe that for moderately stiff chains (), the ratio as expected for ideal chains. For fully flexible chains (), results of the mean square internal distance show remarkable deviations from the freely rotating chain model describing the behavior of ideal chains, while the deviations are diminished as the stiffness of chains increases. For , is in perfect agreement with FRC up to , while a slight deviation occurs for due to the correlation hole effect. Results of the probability distributions of reduced end-to-end distance , and reduced gyration radius for polymer chains in a melt for various values of and for show the nice data collapse, and are described by universal functions, Eqs. (II) and (II), for ideal chains. A detailed investigation of the standard structure factor for single chains in a melt for , , and is also given. Results of presented in a Kratky-plot show that there exists a significant deviation from the Debye function for Gaussian chains at the intermediate values of