Influence of Cohesive Energy on the Thermodynamic Properties of a

Article pubs.acs.org/Macromolecules

Influence of Cohesive Energy on the Thermodynamic Properties of a Model Glass-Forming Polymer Melt Wen-Sheng Xu,*,† Jack F. Douglas,*,§ and Karl F. Freed*,†,‡ †

James Franck Institute and ‡Department of Chemistry, The University of Chicago, Chicago, Illinois 60637, United States § Materials Science and Engineering Division, National Institute of Standards and Technology, Gaithersburg, Maryland 20899, United States ABSTRACT: Monomer chemical structure and architecture represent the most important characteristics of polymers that affect basic molecular parameters (such as the microscopic cohesive energy parameter ϵ and chain persistence length) and that correspondingly govern the bulk physical properties of polymer materials. Here, we focus on elucidating how the microscopic parameter ϵ influences the bulk thermodynamic properties of polymer melts by using molecular dynamics simulations for a standard coarsegrained bead−spring model of unentangled polymer melts under both constant volume and constant pressure conditions. Basic dimensionless thermodynamic properties, such as the cohesive energy density, thermal expansion coefficient, isothermal compressibility, and surface tension, are found to be universal functions of the temperature scaled by ϵ, and thermodynamic signatures for the onset and end of glass formation are identified based on observable features from the static structure factor. We also find that general trends in the thermodynamics and the characteristic temperatures of glass formation determined from our simulations qualitatively accord with the predictions of the generalized entropy theory of polymer glass formation.

1. INTRODUCTION While a complete microscopic molecular understanding of how monomer chemical and molecular structures control polymer properties remains a challenge,1−3 the alteration of monomer chemical and geometric structure has been long recognized as being capable of causing dramatic changes in the nature of polymer glass formation through their influence on basic molecular parameters, such as the microscopic cohesive energy parameter ϵ and chain persistence length.4,5 Naturally, a number of theoretical4−6 and numerical7−11 studies have addressed the influence of these microscopic molecular parameters on polymer glass formation. The fundamental challenge is the development of a fully general theory that predicts how these molecular parameters influence the nature of polymer glass formation. The generalized entropy theory (GET) of polymer glass formation6 provides a convenient vehicle for studying changes in the character of polymer glass formation induced by alterations of basic molecular parameters. The predictive power of the GET is achieved by merging the lattice cluster theory12,13 for the thermodynamics of semiflexible polymers with the Adam−Gibbs (AG) theory,14 which invokes a relationship between the structural relaxation time and the configurational entropy. The lattice cluster theory employs an extended lattice model with structured monomers and describes the influence of various molecular parameters on polymer thermodynamic properties. The resultant GET thus enables estimates of properties associated with glass formation © XXXX American Chemical Society

(such as the structural relaxation time, characteristic temperatures and fragility of glass formation, and high temperature activation energy) as a function of all thermodynamically relevant molecular parameters and thermodynamic conditions such as pressure, volume, dimension, etc. In particular, GET calculations have indicated that both ϵ and the bending rigidity parameter Eb controlling chain stiffness greatly impact upon the dynamics of polymeric glass-forming (GF) liquids through the influence of these parameters on the configurational entropy.15−19 While the GET explains some general trends observed in GF polymers, this theory offers no microscopic mechanism of polymer glass formation due to its largely thermodynamic nature. Computer simulations provide an ideal tool for extraction of useful molecular level information about the dynamics of polymer glass formation, aspects that are normally difficult to isolate separately in experiment since the modification of monomer structure and chemistry in real polymer materials generally leads to changes in more than one molecular parameter. Moreover, computer simulations are necessary to test certain assumptions in the GET, such as the neglect of the entropy of activation, as suggested by AG.14 Simulation also provides invaluable insight into the molecular dependence of the activation energy at high temperatures Received: July 12, 2016 Revised: October 6, 2016

A

DOI: 10.1021/acs.macromol.6b01503 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

2.1. Model. The present work employs a well-studied model of GF “bead−spring” polymer melts. The system is composed of 200 linear chains, each with 16 beads, a chain length that is well below the entanglement length, but long enough to be considered a polymer.26,27 The nonbonded interactions are modeled by the standard truncated-and-shifted LJ potential28

where relaxation is Arrhenius, a quantity that AG assumed could be determined experimentally. The present and subsequent papers focus on the role of the “cohesive energy parameter” ϵ in determining the thermodynamics and dynamics of a model GF polymer melt, a factor neglected in most previous simulations.10 The reason for this neglect is simple; it is generally presumed that changing the temperature is equivalent to changing ϵ, an assumption readily deduced from models of simple fluids, such as one-component Lennard-Jones (LJ) liquids. However, polymers are generally characterized by multiple interparticle potentials governing their bond potentials, chain rigidity, and intermolecular interactions between the chain segments. Hence, it is not obvious that ϵ can simply be “scaled out” of the problem since all these interactions are probably not influenced by ϵ in the same way. The present paper is restricted to the thermodynamics of polymer melts as a reference point for understanding the dynamics of polymer melts, while the subsequent paper20 emphasizes changes in the collective molecular motion that accompany observed changes in thermodynamic properties. We perform extensive molecular dynamics (MD) simulations to analyze the influence of ϵ on basic thermodynamic and dynamic properties of a coarse-grained model of unentangled polymer melts. As an initial step, the present paper considers polymer melts composed of fully flexible chains. Our simulations compare and contrast both isochoric and isobaric thermodynamic paths to polymer glass formation, a point that has received only limited attention in previous computational work.21 We find that basic dimensionless thermodynamic properties (such as the cohesive energy density, number density, thermal expansion coefficient, isothermal compressibility, and surface tension) become unique functions of the temperature scaled by ϵ. An analysis of the static structure factor shows that the onset condition for glass formation is identified by a Hansen−Verlet “freezing” criterion22 that holds for the entire range of ϵ investigated. By exploring the “hyperuniformity parameter”, defined as the ratio of the value of the static structure S(q) at zero wave vector to the height of the first peak of S(q),23−25 we demonstrate that our model polymeric GF fluid progressively approaches a hyperuniform solid state upon cooling. As discussed below, hyperuniformity quantifies the extent to which particles are “jammed” in a global sense, and this concept has recently attracted much interest in characterization of high density packings of hard particles.23−25 Importantly, we show that the temperature demarking the end of glass formation may be estimated from the hyperuniformity parameter. We thus obtain the thermodynamic criteria governing the temperature range of glass formation. Motivated by the GET,6 we additionally investigate the variations of the characteristic temperatures of glass formation with ϵ. Our simulations indicate that all characteristic temperatures increase with ϵ, in good agreement with the GET. Finally, we explore the broadness of glass formation by analyzing the ratios of characteristic temperatures. The broadness of glass formation is nearly independent of ϵ but becomes diminished when shifting the thermodynamic path from isochoric to isobaric, implying that the thermodynamic path strongly affects the nature of polymer glass formation.

12 ⎧ ⎪ 4ϵε[(σ / r ) − (σ /r )6 ] + C(rcut), r < rcut ULJ(r ) = ⎨ ⎪ r ≥ rcut ⎩ 0,

(1)

where r is the distance between the centers of two beads, σ is the effective diameter of the beads, ε sets the energy scale of the system, and the constant C(rcut) ensures that ULJ goes smoothly to zero at the cutoff distance rcut, which we take to equal 2.5σ. Note that the parameter ϵ in eq 1 controls the strength of the nonbonded attractive interactions between the beads, and a similar parameter is an essential feature of the polymer model in the GET of polymer glass formation,6,15−19 thereby motivating the present investigation of the influence of ϵ on polymer glass formation. This ϵ is directly related to the cohesive interaction strength, internal pressure, and heat of vaporization in liquids. In particular, we show in subsection 3.1 that the cohesive energy density, calculated directly from simulations, increases with ϵ for fixed temperatures, so that ϵ can reasonably be termed the “cohesive energy parameter” as its counterpart in the GET.6 The maintenance of bond connectivity between neighboring beads in the chain is ensured by using a harmonic spring potential Uharm(r ) =

1 k b(r − r0)2 2

(2)

where the parameters are taken as kb = 2000ε/σ and r0 = σ so that crystallization is inhibited and the formation of an amorphous material is facilitated.9,29,30 We also explored glass formation at constant pressure in a similar model of polymer melts, where bond connectivity is instead maintained by the finitely extensible nonlinear elastic (FENE) potential.26,27 The results for the FENE bond potential are similar to those for the harmonic bond potential, so the present paper only presents the results for the harmonic bond potential. Some simulations include also a torsional potential,31,32 and moreover, real polymers often contain side groups possessing different degrees of chain stiffness. Hence, our polymer model is clearly rather idealized, but it nonetheless provides a minimal model baseline for investigating glass formation in polymers. The GET incorporates all aforementioned factors, and calculations using the GET indicate that the incorporation of side groups is very important for recovering trends observed in glass formation of real polymers.6 The present work starts from the minimal coarse-grained model described above; additional interactions will be included sequentially to gauge their individual importance. Physical quantities are reported in standard reduced LJ units. All beads have the same mass m, and length and time are given 2

in units of σ and mσ 2/ε , respectively. These reduced units can be transformed to physical units relevant to laboratory measurements on real polymer materials. If we take the size of a chain segment to have the typical magnitude of about 1−2 nm, measure time in picoseconds, and take the energy scale ε of nonbonded bead interactions to be on the order of 1 kJ/mol,

2. MODEL AND SIMULATION DETAILS This section describes the model of polymer melts used in our simulations along with details about molecular dynamics simulations. B


Article

Macromolecules our model would correspond to a typical polymer such as polystyrene whose glass transition temperature for a moderate chain length is Tg ≈ 100 °C. 2.2. Simulation Details. Extensive MD simulations have been performed in three dimensions using the HOOMD-blue simulation package.33−35 Periodic boundary conditions are employed in all directions. The influence of cohesive interaction strength is explored by considering a large range of ϵ. For each ϵ, the present work investigates a wide range of temperatures (T) covering both the high T Arrhenius and low T super-Arrhenius regimes. We also explore the influence of the thermodynamic path on the glass formation of polymer melts for each ϵ by performing simulations under both constant volume (V) and constant pressure (P) conditions. Simulations for each T are performed along an isochoric path at a fixed number density of ρ = N/V = 1 in the NVT ensemble, where N = 3200 is the total number of beads. The pressure is fixed at P = 0 along the isobaric path. This process begins with simulations in the NPT ensemble to determine the desired average density for a given T, and then equilibration and production runs are performed at the target density in the NVT ensemble. This procedure ensures that the average pressure of the system remains approximately constant for all T along the isobaric path. The temperature in the NVT simulations is maintained by the Nosé−Hoover thermostat,36,37 which is implemented in HOOMD-blue via the Martyna−Tobias−Klein equations of motion,38 and the Martyna−Tobias−Klein barostat−thermostat38 is employed in the NPT simulations, where the integration is performed in a cubic box at constant hydrostatic pressure by simultaneously rescaling all three lengths of the simulation box. A time step of Δt = 0.002 is used in both the NVT and NPT simulations. The polymer fluid in all simulations is fully equilibrated before the production runs to strictly avoid nonequilibrium phenomena associated with glass formation. The simulation time for equilibration is generally over 10−100 times larger than the structural relaxation time τα (see the definition of τα in subsection 3.4). Equilibration is also confirmed by monitoring certain thermodynamic quantities (such as potential energy and pressure) during the production runs. Moreover, four independent runs are performed for each state point to obtain reliable results and improve the statistics. Uncertainties shown in certain figures represent standard deviations over four independent runs.

The solubility parameter δ is a fundamental parameter in many engineering applications, and the interaction parameter χ for liquid mixtures is often estimated from measured values of ΠCED for each component.40 From the standpoint of the dynamics of fluids, ΠCED, and closely related properties described below, have been argued to be related to the high temperature activation energy for diffusion and shear viscosity in fluids,41−45 a relation that has been supported by a substantial body of evidence. Moreover, ΠCED has been found to highly correlated with the glass transition temperature Tg of polymer melts from considerations of many materials.46−48 Unfortunately, there is no rigorous development of transition state theory appropriate to describe the dependence of activation free energy parameters on the molecular parameters of polymer fluids, even in the high temperature fluid regime where simple activated transport clearly describes the liquid dynamics phenomenologically. Recent work has addressed this problem by calculating the activation free energy parameters for polymers as a function of chain length based on atomistic modeling,49 but much basic work is necessary to rigorously establish a transition state theory of liquids as a foundation for a fundamental theory of the dynamics of GF liquids. Given the evident importance of ΠCED from both theoretical and polymer engineering perspectives, we first establish the relationship between our molecular parameter ϵ and ΠCED to demonstrate that this quantity provides a measure of the “strength” of molecular cohesion. The magnitude ΠCED for many materials can be tuned over a wide range through the presence of charged and dipole species50−54 or through the addition of salts and/or nanoparticle additives. We expect our calculations for variable ϵ to offer insight into diverse polymeric materials, despite the complexity of the intermolecular interactions involved in general. Our specific model calculations, however, are restricted to polymers interacting with relatively short-ranged van der Waals interactions so that some caution should be made for systems with large ΠCED where long-range interactions are present. The cohesive energy density can be estimated from the intermolecular part of the internal energy for polymeric liquids, as noted in ref 39. Therefore, we calculate ΠCED from our simulations according to the definition

3. RESULTS AND DISCUSSION This section begins by exploring the cohesive energy density of polymer melts, followed by a discussion of the thermodynamic properties that are more relevant to glass formation, such as the density, thermal expansion coefficient, isothermal compressibility, and surface tension. This section then investigates the variations of characteristic temperatures of glass formation with the cohesive energy parameter and measures of the broadness of glass formation. 3.1. Cohesive Energy Density. The cohesive energy density (ΠCED) provides a basic measure of the strength of attractive intermolecular interactions in liquids. Accordingly, ΠCED and related thermodynamic properties have central significance for both the thermodynamics and dynamics of polymer melts. From a thermodynamic standpoint, ΠCED is related to the solubility parameter δ via the relation39

where Ψinter is the total intermolecular potential energy and V is the volume of the simulation box. Since the present work is not intended to investigate polymer melts interacting specifically with a LJ potential, Ψinter and other potential energies are not “corrected” to account for the potential cutoff. The truncatedand-shifted LJ potential is indeed common for modeling polymer melts in MD simulations.9,29,30 Note that ΠCED is related to the internal pressure and the heat of vaporization of volatile fluids and that internal pressure is central to the modeling of the activation energy by Allen and co-workers.55 Figure 1 shows how the computed ΠCED varies with ϵ at a fixed temperature of T = 2 at constant V and constant P. We see that ΠCED grows nearly linearly with ϵ at constant V, although this scaling becomes somewhat nonlinear in ϵ under constant P conditions probably because the volume varies with ϵ along the isobaric path. As an illustration, the dotted line in Figure 1 shows that the variation of ΠCED with ϵ at constant P can be well described by an empirical quadratic function, ΠCED = b0 + b1ϵ + b2ϵ2, where the fitting parameters are provided in

ΠCED = δ 2

ΠCED = −

(3) C

Ψinter V

(4)


Article

Macromolecules

Sanchez56 has recently emphasized that universal scaling relationships for the thermodynamic properties in liquids require the introduction of appropriate dimensionless variables, which are not necessarily the critical parameters corresponding to the liquid−vapor transition as in the van der Waals equation of state. We emphasize that the presence of universal relationships between the reduced thermodynamic properties and the reduced temperature T/ϵ do not imply a trivial influence of ϵ on polymer glass formation since dynamic properties do not scale with T/ϵ in a similar fashion, as discussed in the following paper.20 For completeness, Appendix A discusses the temperature dependence of the nonbonded potential energy Ψnb in our model polymer melt, a quantity related to both the internal energy and specific heat. A universal relationship likewise appears between Ψnb/ϵ and T/ϵ. While the specific heat is one of the most commonly measured thermodynamic properties of polymer melts, this property is one of the least understood properties from a theoretical standpoint. The presence of multiple molecular contributions to the specific heat, reflecting the collective modes and the complex nature of the intermolecular potentials in the fluids, makes the interpretation of the specific heat from simulations difficult. Consequently, computations of the specific heat by MD simulation are very limited for polymers. Since the issues involved are quite subtle, technical, and ultimately unrewarding with respect to the prospect of comparing to experimental measurements for the specific heat, we discuss these issues in Appendix A. 3.2. Equation of State and Thermal Expansion Coefficient. The density as a function of pressure and temperature, or the “equation of state”, represents another fundamental quantity for characterizing the thermodynamics of liquids. Figure 3a displays the temperature variation of the number density ρ for various ϵ at a constant pressure of P = 0. As expected, ρ increases upon cooling for fixed ϵ, and elevating ϵ leads to a growth in ρ for fixed T (see also the inset to Figure 3b). These trends are fully consistent with physical intuition and the previous prediction from the GET.17 In the low T regime, ρ increases nearly linearly with T and extrapolates to ρ* as T → 0, i.e., ρ = ρ* − c(T/ϵ), where c is a fitting parameter. (Sanchez has recently offered an interesting interpretation of ρ*.57) Figure 3b displays ρ/ρ* as a function of T/ϵ, where all data for ρ collapse onto a master curve. The GET likewise predicts that the polymer filling fraction becomes a unique function of T/ϵ for model polymer melts over the entire liquid state regime (data not shown). Sanchez has also emphasized the striking linearity of ρ in T for nonvolatile liquids and the physical meaning of ρ* in his recent reduced variable description of the thermodynamics of liquids.56 The linear T dependence of ρ implies that T and ρ can be taken as essentially interchangeable thermodynamic variables, a property of some importance in comparing theoretical models of the dynamics of GF liquids to measurements. Specifically, when ρ is linear in T, the well-known Vogel−Fulcher− Tamman (VFT) relation58−60 for the α-structural relaxation time τα

Figure 1. Cohesive energy density ΠCED as a function of ϵ for a fixed temperature of T = 2 at constant V (circles) and constant P (squares). Dashed and dotted lines indicate the fits of data for ΠCED to ΠCED = a0 + a1ϵ at constant V and ΠCED = b0 + b1ϵ + b2ϵ2 at constant P, where the fitting parameters (a0, a1) and (b0, b1, b2) are determined to be (−1.74, 5.05) and (−3.52, 4.07, 0.35), respectively.

the caption of Figure 1. Despite the modest differences in ΠCED between constant V and constant P conditions, Figure 1 clearly motivates referring to ϵ as a “cohesive energy parameter”. Figure 2 examines the temperature dependence of ΠCED over a range of T that encompasses both the Arrhenius and super-

Figure 2. (a) Cohesive energy density ΠCED as a function of T for various ϵ [as indicated in (b)] at constant V (dashed lines) and constant P (dotted lines). (b) ΠCED/ϵ as a function of T/ϵ for various ϵ at constant V (dashed lines) and constant P (dotted lines).

Arrhenius regimes of glass formation. As expected, ΠCED increases with decreasing T for each ϵ. In addition, the T dependence of ΠCED becomes much stronger at constant P than at constant V, reflecting the sensitivity of ΠCED to changes in V. Despite the apparently quantitative influence of ϵ on ΠCED, Figure 2b reveals the presence of a universal relationship between the reduced cohesive energy density, ΠCED/ϵ, and the ϵ-scaled temperature, T/ϵ. This universal relationship clearly depends on the thermodynamic path. Note that data for ΠCED do not become a universal function of T/ϵ. Interestingly,

⎛ DT0 ⎞ τα = τ0 exp⎜ ⎟ ⎝ T − T0 ⎠

(5)

becomes “equivalent” to the Doolittle expression61 D


Article

Macromolecules

The thermal expansion coefficient αP is then readily computed from the cubic fits. Motivated by our recent theoretical work based on the GET,64 we consider a dimensionless (or reduced) thermal expansion coefficient TαP, which turns out to be more useful than αP in the description of polymer glass formation. Figure 4 presents TαP as a function of T/ϵ, where a master

Figure 4. Reduced thermal expansion coefficient TαP as a function of T/ϵ for various ϵ at constant P. The dashed line indicates the fit of data for TαP with T/ϵ < 0.8 to the equation TαP = d(T/ϵ), where the fitting parameter is determined to be d = 0.30. The inset depicts TαP as a function of T for various ϵ.

Figure 3. (a) Number density ρ as a function of T for various ϵ at constant P. Solid lines indicate the fits of data for ρ to a cubic polynomial. (b) ρ/ρ* as a function of T/ϵ for various ϵ at constant P. The dashed line indicates the fit of data for ρ with T/ϵ < 0.8 to the equation ρ = ρ* − c(T/ϵ), where the fitting parameters (ρ*, c) are determined to be (1.132, 0.285). The inset to (b) depicts the ϵ dependence of ρ for T = 2. Increasing ϵ tends to elevate the density, similar to increasing the pressure in nonassociating fluids, in the model polymer melt investigated in the present work.

⎛ Av ⎞ τα = τ0 exp⎜ ⎟ ⎝ v − v0 ⎠

curve is apparent, especially at low T. A linear equation, TαP = d(T/ϵ), applies at low T (see the dashed line in Figure 4), where the adjustable parameter is determined to be d = 0.30. Note that a master curve fails to exist between αP and T/ϵ. Simha and Boyer65 have found that the reduced thermal expansion coefficient at the glass transition temperature Tg, i.e., TgαP(Tg), is approximately constant for many polymers. Appendix B demonstrates that Tg/ϵ depends weakly on ϵ in our simulations, suggesting that TgαP(Tg) is nearly constant for various ϵ. By extrapolating TαP to Tg with the linear equation TαP = d(T/ϵ), we find TgαP(Tg) ≈ 0.11. For comparison, Simha and Boyer65 have noted that many polymers conform to the empirical rules, Tg(αP,L − αP,G) = 0.11 and TgαP,L = 0.16, where the subscripts “L” and “G” denote the liquid state above Tg and the nonequilibrium glassy state below Tg, respectively. The present simulation results are limited to the equilibrium liquid regime where αP,L ≈ αP(Tg) so that the second Simha− Boyer rule is more comparable to our estimate from simulations, TgαP(Tg) ≈ 0.11, for the family of coarse-grained models investigated. A similar relationship holds in the context of the melting of crystalline materials, where the melting temperature replaces Tg.66 Recently, we have shown based on the GET64 that TgαP(Tg) varies with fragility, suggesting that this dimensionless expansion coefficient should not be “universal”. 3.3. Static Structure Factor, Isothermal Compressibility, and Surface Tension. Additional basic properties characterizing the thermodynamics of fluids can be obtained from the static structure factor

(6)

where A = D(ρ* − ρ0)/ρ0, ρ0 ≡ ρ(T0), v is the specific volume, and v0 ≡ v(T0). The quantity D in eq 5 is a material specific constant characterizing the “fragility” of GF liquids (discussed below), and T0 is the temperature at which τα extrapolates to infinity, an extrapolation does not imply that τα actually diverges at T0 as the VFT expression is restricted to a temperature regime above Tg. (Estimates of T0 are given as a function of ϵ below.) Equation 6 is a statement regarding a relation between τα and a thermodynamic property, ρ, and does not imply that any general relation between mobility and density applies at a local scale, a common misconception in modeling GF liquids.62 The thermal expansion coefficient may be directly obtained from the density

αP = −

1 ∂ρ ρ ∂T

P

(7)

The empirical Tait equation63 for ρ essentially becomes a quadratic function of T at zero pressure, and this equation has often been shown to quantitatively describe experimental data for the EOS. While a quadratic function is found to capture the T dependence of ρ in our simulations, our simulation data can be described by a cubic function [i.e., ρ = ∑3i=0ciTi with ci (i = 0, ..., 3) being the adjustable parameters] more satisfactorily (see the solid lines in Figure 3a). A linear growth of ρ with T is evident in the low T regime, as discussed earlier. Nevertheless, a cubic function applies in the entire regime of T investigated.

N

S(q) =

N

1 ⟨∑ ∑ exp[−iq·(rj − rk)]⟩ N j=1 k=1

(8)

where i = −1 , q = |q| is the wavenumber, rj is the position of particle j, and ⟨...⟩ denotes the usual thermal average. Illustrative examples are provided in Figure 5 at both constant V and constant P. The static structure factor S(q) describes the E


Article

Macromolecules

Figure 6. Height Sp of the first peak of S(q) at Tl as a function of ϵ at constant V and constant P. Dashed lines indicate averages of each data set over ϵ, and the dotted line indicates the Hansen−Verlet freezing value, Sp = 2.85.

formation,68 another property relating to the dynamics of our model polymer melt, which is discussed below and in the subsequent paper.20 Figure 7a shows that Sp progressively Figure 5. Static structure factor S(q) for various temperatures T for ϵ = 1 at (a) constant V and (b) constant P. Dotted lines indicate linear fits for S(q) at q < 1. Note that the strong increase of S(0) at elevated T under constant P conditions suggests an incipient liquid−vapor transition, where S(0) tends to diverge.

mean correlations in the positions of the segments in the polymer fluid, and the dynamic structure factor describes the rate at which these correlations decay due to spontaneous density fluctuations at equilibrium. The α-structural relaxation time describes the rate at which the relaxation decays as a consequence of the diffusive motion of the polymer segments in the melt. The static structure factor S(q) exhibits mild variations upon cooling in GF liquids, but it enables the determination of at least two important quantities, i.e., the height Sp of the first peak of S(q) and the limiting value S(0) of S(q) as q → 0. The peak height Sp has been utilized as a rough criterion for freezing, namely, the Hansen−Verlet freezing rule,22,67 where freezing occurs when Sp reaches about 2.85. In physical terms, when Sp is large, the particles begin to become “jammed” into their local environment, and a tendency toward local particle localization can be expected near this critical condition Sp ≃ 2.85. However, this localization is transient in equilibrium liquids because the inertia of the polymer segments ultimately allows them to escape from their “cages”. Below, we confirm this interpretation of a thermodynamically defined “onset condition” for glass formation by determining the temperature Tl at which particle localization or “caging” first emerges (see the determination of Tl in subsection 3.4). Figure 6 shows that Sp at the localization temperature Tl is generally quite close to the empirical Hansen−Verlet freezing value (i.e., 2.85) for all ϵ at both constant V and constant P, confirming the hypothesis that T l corresponds to an “amorphous freezing temperature”. This temperature plays an essential role in the development of a universal description of the T dependence of τα, the primary topic of the subsequent paper, which focuses on the dynamics of GF liquids. We note in passing that the T dependence of the position of the first peak has been suggested to correlate with the fragility of glass

Figure 7. (a) Height Sp of the first peak of S(q) and (b) S(0) as a function of T/ϵ for various ϵ at constant V and constant P. Dashes lines in (a) indicate the fits of data for Sp to the equation Sp = e0 + e1(ϵ/T), where the fitting parameters (e0, e1) are determined to be (2.13, 0.73) at constant V and (0.69, 1.39) at constant P, respectively. Dotted lines in (b) indicate the fits of data for S(0) at low T to the power law S(0) = f(T/ϵ)δS, where the fitting parameters (f, δS) are determined to be (0.02, 0.56) at constant V and (0.06, 1.88) at constant P, respectively.

grows upon cooling at both constant V and constant P, consistent with what we might naturally expect with the onset of solidification. As found for the density, the dimensionless property Sp obeys a universal relationship in T/ϵ, where Sp = e0 + e1(ϵ/T) with the fitted values of e0 and e1 indicated in the caption of Figure 7. We next consider a measure of the magnitude of fluctuations in the density, S(0). F


Article

Macromolecules

δh are given in the caption of Figure 8. From this universal curve, the temperature at which h = 6(10−3) is estimated to be about T/ϵ = 0.3 at constant P, a value that is comparable to the characteristic VFT temperature T0 (T0/ϵ ≈ 0.35) determined from τα (see subsection 3.4). We may thus estimate both the “onset” and “end” of the glass-formation process from thermodynamic properties of GF liquids, at least under constant P conditions. On approaching the glassy state, the liquid evidently shows an opposite tendency in comparison to approaching a liquid critical point, in the sense that the compressibility tends to vanish rather than diverge. We emphasize that while the evidence points to a thermodynamic transition underlying glass formation, no phase transition involving the singular changes in the free energy arises as temperature or density is varied. The glass transition is a “rounded” transition, and “fragility” corresponds to the extent of rounding.74 The surface tension γ is another basic thermodynamic property of polymer melts that is crucially important to engineering practice in the area of fiber and film formation and stability, and this property is also a valuable source of thermodynamic information relating to the characterization of polymer melts. Sanchez75 has derived a remarkable relation between γ and the fluid density and compressibility based on continuum Ginzburg−Landau theory

It is well-known that S(0) for an equilibrium fluid is proportional to the isothermal compressibility S(0) = ρkBTκT

(9)

where κT has the thermodynamic definition

κT = −

1 ∂V V ∂P

T

(10)

Since the box length L in simulations is necessarily finite, the smallest wave vector is 2π/L, implying that S(0) must be obtained by extrapolating S(q) to the zero wave vector. Following ref 69, S(0) is estimated by fitting the data for S(q) at small q (i.e., q < 1) to a linear equation, as exhibited by dotted lines in Figure 5. Equation 9 indicates that S(0) is the isothermal compressibility divided by its ideal gas value. This reduced compressibility S(0) is a central quantity in a model of glass formation developed by Schweizer and co-workers,70−72 and S(0) has been indicated to be of great importance in predicting the fragility of glass formation within the GET.64 Figure 7b demonstrates that S(0) at constant V can be described by a power law in T/ϵ, S(0) = f(T/ϵ)δS, where the fitted parameters f and δS are specified in the caption of Figure 7. This scaling only becomes established in the low T regime at constant P, although a maser curve between S(0) and T/ϵ remains in the high T regime. Notice that the variation of S(0) with T is stronger at constant P than at constant V. Torquato and co-workers23−25 have defined the ratio h of S(0) to the peak height Sp to be a measure of fluid “hyperuniformity”, i.e., h = S(0)/Sp, a measure of the extent to which density fluctuations are suppressed at large scales. The hyperuniformity parameter h then provides a quantitative measure of the degree to which the material is globally “jammed”. In an ideal hyperuniform material, h = 0, but a material is “effectively hyperuniform”73 when h is on the order of 10−3−10−4. Figure 8 presents the hyperuniformity parameter

⎛ ρ ⎞1/2 γ = A 0⎜ ⎟ ⎝ κT ⎠

(11)

where A0 is a nearly universal constant for small molecule and polymer liquids.75,76 This relation is supported by observations for numerous fluids, and thus, this relation in conjunction with the relations for ρ and κT allows us to estimate a universal relation for γ as a function of the reduced temperature. Combining eq 11 and our simulation results noted above for ρ and S(0) implies that the surface tension obeys the relation γ 1/2

A 0(kBT )

=

ρ S(0)1/2

(12)

We then expect the reduced surface tension to obey its own universal relation with the reduced temperature, as evidenced in Figure 9. This interesting relation requires conformation by both experiment and simulation. 3.4. Characteristic Temperatures of Glass Formation. As noted before, the GET of glass formation,6 an extension of the AG model14 to include a precise calculation of the configurational entropy rather than the use of experimental estimates of the configurational entropy, describes the glass transition as the consequence of a broad underlying thermodynamic transition exhibiting multiple characteristic temperatures, as observed for rounded transitions.74 These temperatures include the “onset temperature” Tl, a “crossover temperature” Tc, denoting the middle of the transition where polymer melts change from a viscous to rubbery consistency,77 and the “end temperature” T0 of the glass transition, where the structural relaxation time extrapolates to infinity. We again emphasize that no true divergence is expected or implied at T0 by the VFT expression or the AG/GET models6,14 since the measurements are limited to a temperature range above Tg where equilibrium measurements are possible. We next estimate these characteristic temperatures from our simulation data following established methods.21,62,78−82

Figure 8. Hyperuniformity parameter h as a function of T/ϵ for various ϵ at constant V and constant P. Dash-dotted lines indicate the fits of data for h at low T to the power law h = g(T/ϵ)δh, where the fitting parameters (g, δh) are determined to be (0.01, 0.89) at constant V and (0.03, 2.69) at constant P, respectively.

h as a function of T in our model polymer melt at both constant V and constant P. A progressive approach to the hyperuniform state appears upon cooling. For both constant V and P, the dimensionless property h can be described by a master curve in T/ϵ, and in particular, h can be described by a power law h = g(T/ϵ)δh in the low T regime, where the fitted parameters g and G


Article

Macromolecules

estimation of Tl by the method indicated above is evidently difficult by experiment, but the condition that Sp approaches a characteristic value near the Hansen−Verlet freezing value 2.85 seems to offer a simple method of estimating Tl experimentally, as discussed in subsection 3.3. This method requires further validation for other GF liquids before this relation is applied uncritically, however. While the critical value of Sp at the onset temperature is expected to depend on particular GF liquids, this type of criterion should be useful within a class of materials. Since the determinations of other characteristic temperatures (such as Tc and T0) of glass formation involve the structural relaxation time τα, we now describe how τα is estimated from our simulations. Conventionally, τα is obtained in simulations first by computing the “dynamic structure factor” and “intermediate scattering function”. In particular, we consider the self-part of the intermediate scattering function

Figure 9. Theoretical prediction of the reduced surface tension, γ/ A0(kBT)1/2 = ρ/S(0)1/2, as a function of T/ϵ for various ϵ at constant V and constant P based on the theory of Sanchez75 and our simulation results.

N

Fs(q , t ) =

First, the onset temperature Tl corresponds to the emergence of caged or localization particle motion. Previous simulations21,62,82 have estimated Tl as the temperature at which particle caging first emerges. Caging is evidenced in the logarithmic derivative [i.e., ∂ ln ⟨r2(t)⟩/∂ ln(t)] of the meansquared displacement (MSD), which begins to develop a minimum around 1 ps at Tl. Previous work21,62,82 termed this characteristic temperature as the Arrhenius temperature TA since the temperature at which the dynamics starts to be nonArrhenius often emerges close to this temperature, but our following paper indicates that these temperatures can differ from each other, so a different designation of the onset temperature is adopted here. We determine the onset temperature Tl for various ϵ on the basis of MSD. Figure 10 provides examples for illustrating this procedure. The

1 ⟨∑ exp{−iq·[rj(t ) − rj(0)]}⟩ N j=1

(13)

where rj(t) is the position of particle j at the time t. The wavenumber is chosen to be q = 7.1 in this work,83 which is close to the first peak of the static structure factor S(q). Figure 11 displays the time dependence of Fs(q, t) for various T at constant V for ϵ = 1. Following previous work,21,62,78−82 the structural relaxation time τα is defined as Fs(q, t = τα) = 0.2.

Figure 11. Self-intermediate scattering function Fs(q, t) as a function of the time t for various temperatures T for ϵ = 1 at constant V. The wavenumber is chosen to be q = 7.1, which is close to the position of the first peak of S(q). Circles indicate the positions of Fs(q, t = τα) = 0.2, a common convention for defining the structural relaxation time τα. Note that t is given in terms of the reduced units defined in subsection 2.1.

The middle of the glass transition is characterized by a “crossover” temperature Tc, at which a number of subtle changes arise in the properties of polymer fluids. For example, Flory and co-workers84,85 have identified a “third-order transition” near this temperature, and many experimental works45,77,86−90 have emphasized that this “crossover temperature” undoubtedly has significance for polymer processing. The characteristic temperature Tc of glass formation is normally estimated experimentally by fitting τα to a power-law equation τα = τc(T − Tc)−γc

(14)

where τc, Tc, and γc are fitting parameters. This power-law scaling is first inspired by the mode-coupling theory of glass formation,91 and Tc is thus often termed the “mode-coupling temperature”. This power-law description of τα accords with many GF liquids but applies only for a limited temperature

Figure 10. Logarithmic derivative ∂ ln ⟨r2(t)⟩/∂ ln(t) of the meansquared displacement as a function of t for various T for ϵ = 1 at (a) constant V and (b) constant P. The curves marking the “localization temperature” Tl, at which ∂ ln ⟨r2(t)⟩/∂ ln(t) starts to display a local minimum at a time on the order of 1 ps, are indicated as dashed lines. H


Article

Macromolecules

Figure 12. Characteristic temperatures as a function of ϵ at constant V and constant P. (a) Localization temperature Tl. (b) Crossover temperature Tc. (c) Glass transition temperature Tg. (d) VFT temperature T0. Lines indicate the fits of data for characteristic temperatures to the equation Tx = Ix + Jxϵ (x = l, c, g, or 0), where the fitting parameters Ix and Jx are summarized in Table 1. It is apparent that ϵ is a fundamental parameter influencing the dynamical parameters of glass formation, a quantity that may be controlled through the choice of polymer chemistry.

Table 1. Summary of the Fitted Parameters Ix and Jx (x = l, c, g, or 0) for the ϵ Dependence of Various Characteristic Temperatures Shown in Figure 12 at Constant V and Constant P condition

(Il, Jl)

(Ic, Jc)

(Ig, Jg)

(I0, J0)

constant V constant P

(0.33, 1.02) (−0.02, 0.71)

(0.06, 0.35) (0.03, 0.40)

(0.04, 0.29) (0.02, 0.36)

(0.03, 0.25) (0.01, 0.34)

regime above the fitted Tc, but below Tl. We note that the fitted Tc probably does not have anything to do with the predictions of mode-coupling theory since the predicted magnitudes are not close to experimental values for systems where these predictions can be checked.92 Moreover, the same power-law form is predicted to occur in the observed T range by the GET,6 along with the values of the scaling exponent γc. Within the GET, Tc is precisely determined by an inflection point in the product of the configurational entropy density sc times T, a well-defined thermodynamic condition. As with all the characteristic temperatures, the fitted results depend somewhat on the range of T chosen. To determine Tc, the power-law fits are performed for a range of τα with 10 ≤ τα ≤ 103 for each ϵ at both constant V and constant P. The exponent γc determined in this way lies between 2.4 and 3.2, which are comparable with those for other model GF liquids91 and with the GET predictions.6 The ϵ dependence of Tc is discussed below and compared to the other characteristic temperatures. Turning to the determination of Tg and T0, data for τα are first fitted to the VFT equation, as mentioned earlier in subsection 3.2. The “fragility” is defined in eq 5 as KVFT ≡ 1/D so that this quantity conveniently increases as the fragility increases. This is a definition that is useful in engineering practice in defining the strength of the temperature dependence of the dynamics in GF liquids. As before, we emphasize that the fitted results are sensitive to the range of T employed, and we fit τα(T) using a similar range in τα for each ϵ. The lower bounds for τα are chosen to be τα(TA), where TA is the temperature at which the T dependence of τα starts to be nonArrhenius, and the upper bounds for τα are just the largest τα

available in our simulations since they are indeed comparable for all ϵ investigated. The glass transition temperature Tg is then estimated from the VFT fits using the common empirical definition of τα(Tg) = 100 s, where the uncertainties are large because of the large extrapolations involved. We map our reduced variables to physical units relevant to real polymer materials by approximating 1 time unit as 1 ps, as explained in subsection 2.1. Figure 12 summarizes the dependence of various characteristic temperatures of glass formation on ϵ at both constant V and constant P. These characteristic temperatures grow linearly with ϵ, in qualitative agreement with the GET predictions for semiflexible polymer melts.15−19 In particular, the ϵ dependence of these characteristic temperatures can be described by the linear relation Tx = Ix + Jxϵ (x = l, c, g, or 0). The fitted values for Ix and Jx are summarized in Table 1. In addition, our simulations indicate that while the characteristic temperatures (including Tc, Tg, and T0) along the isochoric path with ρ = 1 are found to be lower than those along the isobaric path with P = 0 for fixed ϵ, Tl exhibits the opposite trend. Previous calculations based on the GET have indicated that all characteristic temperatures increase with ϵ in a nonlinear fashion and tend to saturate for large ϵ,15−19 a different trend than indicated in the simulation results in Figure 12. This difference probably derives from the fact that the simulations are performed for a polymer model in which the chains are fully flexible. Indeed, our most recent computations within the GET93 have shown that all characteristic temperatures increase with ϵ exactly in a linear fashion for fully flexible chains. This result suggests that simulations for semiflexible chains are I


Article

Macromolecules

trend is consistent with the GET64 and experiments for polymer melts under isochoric and isobaric conditions.94 By comparing the ratios in Table 2 to the predictions from the GET (see Table 1 in ref 6), we find that our model polymer melt belongs to the class of polymer fluids with both flexible backbone and side groups, a result that comes as no surprise since the polymer chains in the present work are modeled as being fully flexible.

necessary to better understand the dependence of characteristic temperatures on ϵ. Inclusion of a bending potential describing the chain rigidity may alter the influence of ϵ on the dynamics of GF polymer melts in both quantitative and qualitative manners because of the altered coupling between the interaction potential governing the cohesive energy and the bending potential controlling the chain rigidity.

Figure 14. Tg/ϵ as a function of ϵ at constant V and constant P.

4. SUMMARY Despite the long-recognized fact that monomer molecular and chemical structures directly determine the molecular parameters (such as cohesive energy and chain stiffness) governing polymer physical properties, a complete microscopic understanding of how these molecular parameters influence polymer glass formation remains elusive. While the GET6 provides a convenient vehicle for systematically studying changes in polymer glass formation induced by alterations of key molecular parameters such as the microscopic cohesive energy parameter, no microscopic information is offered by the theory at the microscopic scale. Computer simulations are particularly useful to extract microscopic information, and surprisingly, the role of cohesive interaction strength in polymer glass formation has been addressed only in a limited way so far in simulations.10 To better understand how the cohesive interaction strength influences polymer glass formation, we perform extensive molecular dynamics simulations for a coarse-grained model of GF polymer melts under both constant volume and constant pressure conditions. The present paper particularly analyzes several basic thermodynamic properties upon cooling. Our simulations indicate universal scaling behavior between dimensionless thermodynamic properties (such as the cohesive energy density, number density, thermal expansion coefficient, isothermal compressibility, and surface tension) and the ϵscaled temperature. An analysis of the static structure factor shows that the onset temperature of glass formation may be estimated by a Hansen−Verlet freezing criterion,22 thereby suggesting a simple method to estimate the onset temperature of glass formation in experiment. We also demonstrate that polymer fluids progressively approach to a hyperuniform glassy state upon cooling. Importantly, we show that the temperature demarking the end of glass formation may be estimated from the magnitude of the hyperuniformity parameter of Torquato and co-workers.23−25 Our work thus indicates thermodynamic signatures corresponding to the onset and end of glass formation. Previous

Figure 13. (a) Nonbonded potential energy per bead Ψnb/N as a function of T for various ϵ [as indicated in (b)] at constant V (dashed lines) and constant P (dotted lines). (b) Ψnb/(Nϵ) as a function of T/ ϵ for various ϵ at constant V (dashed lines) and constant P (dotted lines).

The GET predicts that the ratios of the characteristic temperatures provide a model independent measure of the broadness of glass formation and thus the “fragility” of glass formation in its most fundamental sense of the degree of transition broadening.6 (Engineers also speak of “long” and “short” corresponding to fragile and strong glass formation in relation to the working time with a material at a given cooling rate before it solidifies.) The ratios of the characteristic temperatures are found to be nearly constants (data not shown), meaning that ϵ does not change the fragility within this family of polymer melts. Table 2 summarizes the averaged Table 2. Averaged Ratios of the Characteristic Temperatures over ϵ at Constant V and Constant P condition

Tl/Tc

Tl/Tg

Tl/T0

Tc/Tg

Tc/T0

constant V constant P

3.09 1.69

3.84 1.88

4.36 2.01

1.24 1.11

1.41 1.19

ratios of the characteristic temperature over ϵ at constant V and constant P. The broadness of glass formation clearly becomes reduced upon shifting the thermodynamic path from constant V to constant P conditions, implying that the fragility is larger at constant P than at constant V. We further analyze this feature of glass formation in our following paper, which emphasizes the dynamics of the present model of polymer melts. This general J


Article

Macromolecules experimental works by Flory and co-workers,84,85 and many others,45,77,86−90 have suggested that the crossover temperature Tc of glass formation also has a thermodynamic signature in the temperature dependence of the density, compressibility, or other properties. This characteristic temperature is predicted by the GET to be the inflection point in the product of the configurational entropy and temperature6 but this temperature has been given alternative interpretations in other models of glass formation. Unfortunately, investigations of this low temperature regime of glass formation below Tc by molecular dynamics are quit difficult at present. Nonetheless, we may infer from our simulations that glass formation is consistent with a rounded thermodynamic transition, a type of transition often found in molecular self-assembly,21,74 but there is no evidence for a true phase transition in which the free energy changes singularly with temperature at some critical temperature. Ironically, while the other characteristic temperatures apparently have some thermodynamic “meaning”, the “glass transition temperature” itself is just a useful engineering parameter that is found by locating “kinks” in various quasithermodynamic properties, a feature that implies a small temperature range where it is difficult for a material to achieve equilibrium under some specified rate of cooling. This widely reported parameter correspondingly has no obvious fundamental significance in terms of fluid thermodynamics. We have minimized our discussion of this commonly reported temperature for GF liquids, and we believe that too much emphasis has been given to this temperature whose significance is fraught with difficulties. Finally, we analyze how the cohesive energy parameter ϵ influences the characteristic temperatures of glass formation, motivated by the GET.6 Our simulations indicate that all characteristic temperatures of glass formation increase with ϵ, in good agreement with the GET. Accordingly, an analysis of ratios of the characteristic temperatures permits investigations for the broadness of glass formation. We find that the “broadness” of glass formation is nearly independent of ϵ but becomes reduced when shifting the thermodynamic conditions from a constant volume to a constant pressure, implying that the thermodynamic path has a large influence on the nature of polymer glass formation.

Since the contribution from the bonded interactions to the internal energy depends largely on the type of the bond potential, we mainly consider the total nonbonded potential energy Ψnb and its contribution to the specific heat. Figure 13 displays the temperature dependence of Ψnb/N as decreasing with T at both constant V and constant P. As found for the cohesive energy density, a universal relationship, which depends on the thermodynamic path, appears between Ψnb/(Nϵ) and T/ϵ. The isochoric and isobaric specific heats can be normally determined from the internal energy U and the enthalpy H = U + PV as CV = ∂U/∂T |V and CP = ∂H/∂T |P, respectively. Since our simulations for constant P conditions are performed at P = 0, both CV and CP are simply the first derivative of U with respect to T. (The average kinetic energy is equal to (3/2)kBT from equipartition so this contribution to the specific heat is simply omitted here.) The results in Figure 13b then imply that the nonbonded contribution to the specific heat are unique functions of T/ϵ for various ϵ in the present model of polymer melts. Again, we anticipate that the average potential energy should be sensitive to the intermolecular potentials involved so that the generality of these scalings for other polymer models remains to be investigated. The specific heat of polymer chains, being described rather ideally as chains of volume-excluding and self-attracting harmonic oscillators in a coarse-grained sense, must also exhibit anharmonic interactions from these interactions. It is then not clear to us that the nonbonded contribution to the specific heat can be directly compared to experimental estimates of the specific heat. We really need a more thorough understanding of the vibrational contribution to the specific heat beyond the simple harmonic approximation. Since a positive slope in the T dependence of the specific heat above Tg is typical of polymer materials,99,100 we have a strong suggestion from measurement that the vibrational contribution to the specific heat of polymers might be large. Consistent with this hypothesis, Starr et al.80 found that the configurational entropy of a similar coarsegrained polymer melt to the present model only contributes about 10% of the total entropy (see also ref 98), while the vibrational entropy accounts for the rest. The specific heat of polymers is then primarily vibrational with a configurational “background”. Of course, this is a serious problem in polymeric materials when trying to extract the configurational contribution to the fluid entropy from a specific heat measurement since there is no reason to believe that the anharmonic contributions to polymer materials should be the same as in the low temperature solid state, be it either a crystalline or glassy material. We contrast this situation with small molecule liquids, such as water, where the configurational entropy makes a relatively large contribution to the total fluid entropy.101 It is then no wonder that estimates of the “thermodynamic fragility” from specific heat measurements are more informative about the dynamics of small molecule liquids than polymer liquids.102,103 We can also understand why the excess entropy of polymer materials can extrapolate to zero at a temperature that is quite different from T0. Ironically, the specific heat is one of the most intensively studied properties for polymer materials. Wunderlich and coworkers,99,100 in particular, have made great efforts in tabulating and interpreting data for the specific heat taken for diverse polymers over a wide range of temperatures. This combined measurement and modeling effort has suggested a possible origin of the “anomalous” linear increase in the specific heat for

■

APPENDIX A. ISSUES WITH ESTIMATING THE SPECIFIC HEAT Although the internal energy and specific heat are widely considered thermodynamic properties from an experimental standpoint, these properties are the least understood and studied thermodynamic properties of fluids from a theoretical standpoint.95,96 There is also a surprising dearth of molecular dynamics studies of these properties. The conceptual problem here is that it remains unclear what kinds of normal modes contribute to these properties given the fact that the disorder and anharmonic effects must be prevalent in the fluid state.95 We approach this notoriously difficult problem by simulation estimates of the internal energy and then calculate the specific heat from this property. Simulations also allow to assess the different contributions to the specific heat that arise from bonded and nonbonded interactions and to separate the vibrational contributions from the configurational contributions.97,98 Therefore, the relative importance of these contributions can be examined in polymeric and small molecule liquids. K


Macromolecules

■

many polymers at temperatures above Tg. In particular, the Tarasov model104 indicates that the specific heat has a large vibrational contribution that arises from chain connectivity which is modeled by treating polymers as one-dimensional crystals and the invocation of the Debye theory of the specific heat in one dimension. This simple model has been found to fit data for the specific heat rather well, where the Debye temperature (θ) is found to be rather large105,106 (e.g., 519 K for polyethylene, 713 K for polypropylene, and 1721 K for diamond) because of the large force constants of the chemical bonds of the polymer chains. The specific heat for a crystal scales as CV ∼ T for temperatures below about θ/3, so we have an evident physical rationale for the linear T scaling of the specific heat over a large range of T above Tg. Further work is needed to estimate the validity of the Tarasov model and to calculate the Debye temperature directly from knowledge of the bond vibrational frequencies of polymers.

■

APPENDIX B. GLASS TRANSITION TEMPERATURE

■

AUTHOR INFORMATION

REFERENCES

(1) Ngai, K. L.; Roland, C. M. Chemical Structure and Intermolecular Cooperativity: Dielectric Relaxation Results. Macromolecules 1993, 26, 6824−6830. (2) Kunal, K.; Robertson, C. G.; Pawlus, S.; Hahn, S. F.; Sokolov, A. P. Role of Chemical Structure in Fragility of Polymers: A Qualitative Picture. Macromolecules 2008, 41, 7232−7238. (3) Agapov, A. L.; Wang, Y.; Kunal, K.; Robertson, C. G.; Sokolov, A. P. Effect of Polar Interactions on Polymer Dynamics. Macromolecules 2012, 45, 8430−8437. (4) Gibbs, J. H.; DiMarzio, E. A. Nature of the Glass Transition and the Glassy State. J. Chem. Phys. 1958, 28, 373−383. (5) DiMarzio, E. A.; Gibbs, J. H. Chain Stiffness and the Lattice Theory of Polymer Phases. J. Chem. Phys. 1958, 28, 807−813. (6) Dudowicz, J.; Freed, K. F.; Douglas, J. F. Generalized Entropy Theory of Polymer Glass Formation. Adv. Chem. Phys. 2008, 137, 125−222. (7) Barrat, J.-L.; Baschnagel, J.; Lyulin, A. Molecular dynamics simulations of glassy polymers. Soft Matter 2010, 6, 3430−3446. (8) Kumar, R.; Goswami, M.; Sumpter, B. G.; Novikov, V. N.; Sokolov, A. P. Effects of backbone rigidity on the local structure and dynamics in polymer melts and glasses. Phys. Chem. Chem. Phys. 2013, 15, 4604−4609. (9) Shavit, A.; Riggleman, R. A. Influence of Backbone Rigidity on Nanoscale Confinement Effects in Model Glass-Forming Polymers. Macromolecules 2013, 46, 5044−5052. (10) Xia, W.; Keten, S. Coupled Effects of Substrate Adhesion and Intermolecular Forces on Polymer Thin Film Glass-Transition Behavior. Langmuir 2013, 29, 12730−12736. (11) Xie, S.-J.; Qian, H.-J.; Lu, Z.-Y. The glass transition of polymers with different side-chain stiffness confined in free-standing thin films. J. Chem. Phys. 2015, 142, 074902. (12) Foreman, K. W.; Freed, K. F. Lattice Cluster Theory of Multicomponent Polymer Systems: Chain Semiflexibility and Specific Interactions. Adv. Chem. Phys. 1998, 103, 335−390. (13) Xu, W.-S.; Freed, K. F. Lattice cluster theory for polymer melts with specific interactions. J. Chem. Phys. 2014, 141, 044909. (14) Adam, G.; Gibbs, J. H. On the Temperature Dependence of Cooperative Relaxation Properties in Glass-Forming Liquids. J. Chem. Phys. 1965, 43, 139−146. (15) Stukalin, E. B.; Douglas, J. F.; Freed, K. F. Application of the entropy theory of glass formation to poly(α-olefins). J. Chem. Phys. 2009, 131, 114905. (16) Xu, W.-S.; Freed, K. F. Thermodynamic scaling of dynamics in polymer melts: Predictions from the generalized entropy theory. J. Chem. Phys. 2013, 138, 234501. (17) Xu, W.-S.; Freed, K. F. Influence of Cohesive Energy and Chain Stiffness on Polymer Glass Formation. Macromolecules 2014, 47, 6990−6997. (18) Xu, W.-S.; Freed, K. F. Generalized Entropy Theory of Glass Formation in Polymer Melts with Specific Interactions. Macromolecules 2015, 48, 2333−2343. (19) Dudowicz, J.; Douglas, J. F.; Freed, K. F. Advances in the generalized entropy theory of glass-formation in polymer melts. J. Chem. Phys. 2014, 141, 234903. (20) Xu, W.-S.; Douglas, J. F.; Freed, K. F. Influence of Cohesive Energy on Relaxation in a Model Glass-Forming Polymer Melt. Macromolecules 2016, DOI: 10.1021/acs.macromol.6b01504. (21) Pazmiño Betancourt, B. A.; Douglas, J. F.; Starr, F. W. String model for the dynamics of glass-forming liquids. J. Chem. Phys. 2014, 140, 204509. (22) Hansen, J.-P.; Verlet, L. Phase Transitions of the Lennard-Jones System. Phys. Rev. 1969, 184, 151−161. (23) Torquato, S.; Stillinger, F. H. Local density fluctuations, hyperuniformity, and order metrics. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2003, 68, 041113. (24) Zachary, C. E.; Torquato, S. Hyperuniformity in point patterns and two-phase random heterogeneous media. J. Stat. Mech.: Theory Exp. 2009, 2009, P12015.

In the main text, we demonstrate that the reduced thermal expansion coefficient TαP is a universal function of Tg/ϵ for our simulation model. Figure 14 displays Tg/ϵ as depending very weakly on ϵ at both constant V and constant P. Notably, Tg/ϵ appears to slightly decrease with increasing ϵ. However, given the fact that Tg is estimated by extrapolating the data at much higher T, Figure 14 suggests that TgαP(Tg) stays roughly constant for various ϵ, consistent with the analysis of Simha and Boyer.65 By taking Tg/ϵ ≈ 0.37 at constant P, we find TgαP(Tg) ≈ 0.11 using the linear fit TαP = 0.30(T/ϵ) discussed in the main text. Note that Tg/ϵ at constant V are quite different from those at constant P.

Corresponding Authors

*E-mail: [email protected] (K.F.F.). *E-mail: [email protected] (W.-S.X.). *E-mail: [email protected] (J.F.D.). Present Address

W.-S.X.: Center for Nanophase Materials Sciences, Oak Ridge National Laboratory, Oak Ridge, TN 37831. Notes

The authors declare no competing financial interest.

■

Article

ACKNOWLEDGMENTS

We thank Prof. Salvatore Torquato (Princeton University) for helpful discussions on hyperuniformity and valuable comments on the manuscript, Prof. Francis Starr (Wesleyan University) for helpful conversations, and Dr. Alexandros Chremos (NIST) for useful comments on the manuscript. W.-S.X. is grateful to Prof. Juan J. de Pablo and his group members for providing the opportunity to attend their group meeting while working at the University of Chicago, from which the present work has greatly benefited. We are grateful for the support of the University of Chicago Research Computing Center for assistance with the simulations carried out in this work. This work is supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering under Award DE-SC0008631. L


Article

Macromolecules (25) Atkinson, S.; Zhang, G.; Hopkins, A. B.; Torquato, S. Critical slowing down and hyperuniformity on approach to jamming. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2016, 94, 012902. (26) Kremer, K.; Grest, G. S. Dynamics of entangled linear polymer melts: A molecular-dynamics simulation. J. Chem. Phys. 1990, 92, 5057−5086. (27) Grest, G. S.; Kremer, K. Molecular dynamics simulation for polymers in the presence of a heat bath. Phys. Rev. A: At., Mol., Opt. Phys. 1986, 33, 3628−3631. (28) Frenkel, D.; Smit, B. Understanding Molecular Simulation: from Algorithms to Applications; Academic Press: London, 2002. (29) Riggleman, R. A.; Yoshimoto, K.; Douglas, J. F.; de Pablo, J. J. Influence of Confinement on the Fragility of Antiplasticized and Pure Polymer Films. Phys. Rev. Lett. 2006, 97, 045502. (30) Riggleman, R. A.; Douglas, J. F.; de Pablo, J. J. Tuning polymer melt fragility with antiplasticizer additives. J. Chem. Phys. 2007, 126, 234903. (31) Bulacu, M.; van der Giessen, E. Effect of bending and torsion rigidity on self-diffusion in polymer melts: A molecular-dynamics study. J. Chem. Phys. 2005, 123, 114901. (32) Boland, E. K.; Liu, J.; Maranas, J. K. A molecular picture of motion in polyolefins. J. Chem. Phys. 2010, 132, 144901. (33) Anderson, J. A.; Lorenz, C. D.; Travesset, A. General purpose molecular dynamics simulations fully implemented on graphics processing units. J. Comput. Phys. 2008, 227, 5342−5359. (34) Glaser, J.; Nguyen, T. D.; Anderson, J. A.; Lui, P.; Spiga, F.; Millan, J. A.; Morse, D. C.; Glotzer, S. C. Strong scaling of generalpurpose molecular dynamics simulations on GPUs. Comput. Phys. Commun. 2015, 192, 97−107. (35) HOOMD-blue Web page: http://glotzerlab.engin.umich.edu/ hoomd-blue. (36) Nosé, S. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys. 1984, 81, 511−519. (37) Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. (38) Martyna, G. J.; Tobias, D. J.; Klein, M. L. Constant pressure molecular dynamics algorithms. J. Chem. Phys. 1994, 101, 4177−4189. (39) Maranas, J. K.; Mondello, M.; Grest, G. S.; Kumar, S. K.; Debenedetti, P. G.; Graessley, W. W. Liquid Structure, Thermodynamics, and Mixing Behavior of Saturated Hydrocarbon Polymers. 1. Cohesive Energy Density and Internal Pressure. Macromolecules 1998, 31, 6991−6997. (40) Merk, W.; Lichtenthaler, R. N.; Prausnitz, J. M. Solubilities of fifteen solvents in copolymers of poly(vinyl acetate) and poly(vinyl chloride) from gas-liquid chromatography. Estimation of polymer solubility parameters. J. Phys. Chem. 1980, 84, 1694−1698. (41) Tabor, D. Micromolecular processes in the viscous flow of hydrocarbons. Philos. Mag. A 1988, 57, 217−224. (42) Rennie, A. R.; Tabor, D. Diffusion in Molten Polymers: The Influence of Hydrostatic Pressure. Proc. R. Soc. London, Ser. A 1982, 383, 1−14. (43) Kincaid, J. F.; Eyring, H.; Stearn, A. E. The Theory of Absolute Reaction Rates and its Application to Viscosity and Diffusion in the Liquid State. Chem. Rev. 1941, 28, 301−365. (44) Bershtein, V. A.; Egorov, V. M.; Egorova, L. M.; Ryzhov, V. A. The role of thermal analysis in revealing the common molecular nature of transitions in polymers. Thermochim. Acta 1994, 238, 41−73. (45) Bershtein, V. A.; Ryzhov, V. A. Far Infrared Spectroscopy of Polymers. Adv. Polym. Sci. 1994, 114, 43−121. (46) He, T.-B. On the chain cross-sectional area and motion of macromolecular chains. J. Appl. Polym. Sci. 1985, 30, 4319−4324. (47) Kreibch, U. T.; Batzer, H. Influence of the segment structure and crosslinking on the glass transition Tg. Possibilities of predicting Tg using the values of cohesive energy Ecoh. Angew. Makromol. Chem. 1979, 83, 57−112. (48) Boyer, R. F. Dependence of Tg(K) on the product of the cohesive energy density (CED) and chain stiffness parameter C.infin. Macromolecules 1992, 25, 5326−5330.

(49) Jeong, C.; Douglas, J. F. Mass dependence of the activation enthalpy and entropy of unentangled linear alkane chains. J. Chem. Phys. 2015, 143, 144905. (50) Eisenberg, A.; Farb, H.; Cool, L. G. Glass Transitions in Ionic Polymers. J. Polym. Sci., Part A-2 1966, 4, 855−868. (51) Eisenberg, A.; Matsura, H.; Yokoyama, T. Glass Transitions in Ionic Polymers: The acrylates. J. Polym. Sci., Part A-2 1971, 9, 2131− 2135. (52) Eisenberg, A. Glass Transitions in Ionic Polymers. Macromolecules 1971, 4, 125−128. (53) Eisenberg, A.; Saito, S. Possible Experimental Equivalence of the Gibbs-DiMarzio and Free-Volume Theories of the Glass Transition. J. Chem. Phys. 1966, 45, 1673−1678. (54) Tsutsui, T.; Tanaka, T. Role of electrostatic forces in the glass transition temperatures of ionic polymers. Polymer 1977, 18, 817−821. (55) Allen, G.; Gee, G.; Lanceley, H. A.; Mangaraj, D. Intermolecular forces and chain flexibility in linear polymers. J. Polym. Sci. 1959, 34, 349−354. (56) Sanchez, I. C. Dimensionless Thermodynamics: A New Paradigm for Liquid State Properties. J. Phys. Chem. B 2014, 118, 9386−9397. (57) Sanchez, I. C.; O’Keefe, S.; Xu, J. F. New Zeno-Like Liquid States. J. Phys. Chem. B 2016, 120, 3705−3712. (58) Vogel, H. The law of the relationship between viscosity of liquids and the temperature. Phys. Z. 1921, 22, 645−646. (59) Fulcher, G. S. Analysis of recent measurements of the viscosity of glasses. J. Am. Ceram. Soc. 1925, 8, 339−355. (60) Tammann, G.; Hesse, W. Die Abhängigkeit der viscosität von der temperatur bie unterkühlten flüssigkeiten. Z. Anorg. Allg. Chem. 1926, 156, 245−257. (61) Doolittle, A. K. Studies in Newtonian Flow. I. The Dependence of the Viscosity of Liquids on Temperature. J. Appl. Phys. 1951, 22, 1031−1035. (62) Hanakata, P. Z.; Pazmiño Betancourt, B. A.; Douglas, J. F.; Starr, F. W. A unifying framework to quantify the effects of substrate interactions, stiffness, and roughness on the dynamics of thin supported polymer films. J. Chem. Phys. 2015, 142, 234907. (63) Floudas, G.; Paluch, M.; Grzybowski, A.; Ngai, K. L. Molecular Dynamics of Glass-Forming Systems: Effects of Pressure, 1st ed.; SpringerVerlag: Berlin, 2010. (64) Xu, W.-S.; Douglas, J. F.; Freed, K. F. Entropy Theory of Polymer Glass-Formation in Variable Spatial Dimension. Adv. Chem. Phys. 2016, 161, 443−497. (65) Simha, R.; Boyer, R. F. On a General Relation Involving the Glass Temperature and Coefficients of Expansion of Polymers. J. Chem. Phys. 1962, 37, 1003−1007. (66) Granato, A. V.; Joncich, D. M.; Khonik, V. A. Melting, thermal expansion, and the Lindemann rule for elemental substances. Appl. Phys. Lett. 2010, 97, 171911. (67) Hoffmann, G. P.; Löwen, H. Freezing and melting criteria in non-equilibrium. J. Phys.: Condens. Matter 2001, 13, 9197−9206. (68) Mauro, N. A.; Blodgett, M.; Johnson, M. L.; Vogt, A. J.; Kelton, K. F. A structural signature of liquid fragility. Nat. Commun. 2014, 5, 4616. (69) Marcotte, É; Stillinger, F. H.; Torquato, S. Nonequilibrium static growing length scales in supercooled liquids on approaching the glass transition. J. Chem. Phys. 2013, 138, 12A508. (70) Schweizer, K. S.; Saltzman, E. J. Entropic barriers, activated hopping, and the glass transition in colloidal suspensions. J. Chem. Phys. 2003, 119, 1181−1196. (71) Schweizer, K. S.; Saltzman, E. J. Theory of dynamic barriers, activated hopping, and the glass transition in polymer melts. J. Chem. Phys. 2004, 121, 1984−2000. (72) Saltzman, E. J.; Schweizer, K. S. Universal scaling, dynamic fragility, segmental relaxation, and vitrification in polymer melts. J. Chem. Phys. 2004, 121, 2001−2009. (73) Torquato, S. Personal communication. M


Article

Macromolecules (74) Douglas, J. F.; Dudowicz, J.; Freed, K. F. Does equilibrium polymerization describe the dynamic heterogeneity of glass-forming liquids? J. Chem. Phys. 2006, 125, 144907. (75) Sanchez, I. C. Liquids: Surface tension, compressibility, and invariants. J. Chem. Phys. 1983, 79, 405−415. (76) Bhatia, Q. S.; Chen, J.-K.; Koberstein, J. T.; Sohn, J. E.; Emerson, J. A. The measurement of polymer surface tension by drop image processing: Application to PDMS and comparison with theory. J. Colloid Interface Sci. 1985, 106, 353−359. (77) Hedvat, S. Molecular interpretation of T||, the rubbery-viscous ‘transitiona temperature of amorphous polymers. Polymer 1981, 22, 774−777. (78) Starr, F. W.; Douglas, J. F. Modifying Fragility and Collective Motion in Polymer Melts with Nanoparticles. Phys. Rev. Lett. 2011, 106, 115702. (79) Pazmiño Betancourt, B. A.; Douglas, J. F.; Starr, F. W. Fragility and cooperative motion in a glass-forming polymer-nanoparticle composite. Soft Matter 2013, 9, 241−254. (80) Starr, F. W.; Douglas, J. F.; Sastry, S. The relationship of dynamical heterogeneity to the Adam-Gibbs and random first-order transition theories of glass formation. J. Chem. Phys. 2013, 138, 12A541. (81) Hanakata, P. Z.; Douglas, J. F.; Starr, F. W. Interfacial mobility scale determines the scale of collective motion and relaxation rate in polymer films. Nat. Commun. 2014, 5, 4163. (82) Pazmiño Betancourt, B. A.; Hanakata, P. Z.; Starr, F. W.; Douglas, J. F. Quantitative relations between cooperative motion, emergent elasticity, and free volume in model glass-forming polymer materials. Proc. Natl. Acad. Sci. U. S. A. 2015, 112, 2966−2971. (83) This choice for q is based on the following consideration. It seems to provide the best correspondence between structural relaxation time data from mechanical measurements and those from the self-intermediate scattering function. It is not really known at present why this choice “works”. We plan to explore the relation between the self-intermediate scattering function and the stress relaxation function in the future to better understand this empirical relationship. We note that the temperature dependence of the structural relaxation time will remain qualitatively unchanged if other values of wave vectors are used as long as the choice does not differ too much from the value corresponding to the first peak. (84) Fox, T. G., Jr.; Flory, P. J. Second-Order Transition Temperatures and Related Properties of Polystyrene. I. Influence of Molecular Weight. J. Appl. Phys. 1950, 21, 581−591. (85) Hö cker, H.; Blake, G. J.; Flory, P. J. Equation-of-state parameters for polystyrene. Trans. Faraday Soc. 1971, 67, 2251−2257. (86) Boyer, R. F. Dynamics and thermodynamics of the liquid state (T < Tg) of amorphous polymers. J. Macromol. Sci., Part B: Phys. 1980, 18, 461−553. (87) Boyer, R. F. T|| from the temperature dependence of carbon-13 NMR correlation times τc. J. Polym. Sci., Part B: Polym. Phys. 1988, 26, 893−910. (88) Lacabanne, C.; Goyaud, P.; Boyer, R. F. Thermally stimulated current (TSC) study of the Tg and T|| transitions in anionic polystyrenes. J. Polym. Sci., Polym. Phys. Ed. 1980, 18, 277−284. (89) Boyer, R. F. Pressure dependence of secondary transitions in amorphous polymers. 1. T|| for polystyrene, poly(vinyl acetate), and polyisobutylene. Macromolecules 1981, 14, 376−385. (90) Boyer, R. F.; Miller, R. L. Correlation of liquid-state compressibility and bulk modulus with cross-sectional area per polymer chain. Macromolecules 1984, 17, 365−369. (91) Götze, W. Complex Dynamics of Glass-Forming Liquids: A ModeCoupling Theory; Oxford University Press: Oxford, 2008. (92) Zhang, H.; Zhong, C.; Douglas, J. F.; Wang, X.; Cao, Q.; Zhang, D.; Jiang, J.-Z. Role of string-like collective atomic motion on diffusion and structural relaxation in glass forming Cu-Zr alloys. J. Chem. Phys. 2015, 142, 164506. (93) Xu, W.-S.; Douglas, J. F.; Freed, K. F. Unpublished.

(94) Huang, D.; Colucci, D. M.; McKenna, G. B. Dynamic fragility in polymers: A comparison in isobaric and isochoric conditions. J. Chem. Phys. 2002, 116, 3925−3934. (95) Brazhkin, V. V.; Trachenko, K. Collective Excitations and Thermodynamics of Disordered State: New Insights into an Old Problem. J. Phys. Chem. B 2014, 118, 11417−11427. (96) Granato, A. V. The specific heat of simple liquids. J. Non-Cryst. Solids 2002, 307−310, 376−386. (97) Goldstein, M. Viscous liquids and the glass transition. V. Sources of the excess specific heat of the liquid. J. Chem. Phys. 1976, 64, 4767−4774. (98) Roe, R.-J.; Tonelli, A. E. Contribution of the Conformational Specific Heat of Polymer Chains to the Specific Heat Difference between Liquid and Glass. Macromolecules 1978, 11, 114−117. (99) Gaur, U.; Wunderlich, B. Heat Capacity and Other Thermodynamic Properties of Linear Macromolecules. V. Polystyrene. J. Phys. Chem. Ref. Data 1982, 11, 313−325. (100) ATHAS data base: http://web.archive.org/web/ 20060829233626/http://web.utk.edu/athas/databank/intro.html. (101) Starr, F. W.; Angell, C. A.; Stanley, H. E. Prediction of entropy and dynamic properties of water below the homogeneous nucleation temperature. Phys. A 2003, 323, 51−66. (102) Huang, D.; McKenna, G. B. New insights into the fragility dilemma in liquids. J. Chem. Phys. 2001, 114, 5621−5630. (103) Martinez, L.-M.; Angell, C. A. A thermodynamic connection to the fragility of glass-forming liquids. Nature 2001, 410, 663−667. (104) Tarasov, V. V.; Yunitskii, G. A. Theory of heat capacity of chain-layer structures. Zh. Fiz. Chim. 1965, 39, 2077−2080. (105) Pyda, M.; Bartkowiak, M.; Wunderlich, B. Computation of Heat Capacities of Solids Using a General Tarasov Equation. J. Therm. Anal. Calorim. 1998, 52, 631−656. (106) Jin, Y.; Wunderlich, B. Heat capacities of paraffins and polyethylene. J. Phys. Chem. 1991, 95, 9000−9007.

N


Influence of Cohesive Energy on the Thermodynamic Properties of a

Recommend Documents