Heterogeneous dynamics in linear polymer melts near glass transition temperature

Equilibrium molecular dynamics simulations of a model linear polymer melt is carried out near its glass transition temperature Tc. We have elucidated the role played by intra and inter chain monomers in formation of secondary structures using radial distribution functions on reduction of temperature towards Tc. Our studies on dynamics using intermediate scattering function Fs(k,t) shows step like relaxation, which is a signature of heterogeneous dynamics of the glassy domains that is widely known in earlier studies of glass transition. Further studies of heterogeneity in the dynamics using four-point correlation function x4 s(t) shows the growth of correlated region in this system as temperature is reduced towards Tc.


Introduction
Glass transition and glassy dynamics has been a subject of intense research in last several decades: as glassy systems are abundant and glass transition is ubiquitous in nature [1,2]. Structural glasses are special class of glassy systems which requires structural rearrangement for relaxation near the glass transition temperature T c . Many systems exhibit structural glass transition, examples include systems of polymers, bio-polymers etc.
An important theoretical approach to understand changes in the microscopic dynamics during the glass transition is via studies of the mode coupling theory(MCT) [3][4][5]. The ideal MCT [6,7] has explained the temperature dependency of α-relaxation time, viscosity, and diffusion coefficients etc., above glass transition temperature T c . One of the prominent feature of glass transition is exponential slowing down of the relaxation process; in experiments the glass transition occurs at around viscosity equal to 10 13 poise(or, relaxation time ≈ 100 sec) [8]. Below the crossover temperatures of supercooled liquids, particle motion become increasingly correlated, resulting in formation of cages by neighboring particles. The relaxation of a particle within a cage results in β relaxation and relaxation of the cages classified as α relaxation [9,10]. The structures of these cages around each particle differ which results in dynamic heterogeneity [11,12]. Computational studies of the glass transition in various glass forming model liquids are successful in explaining many microscopic features such as microscopic relaxation, dynamic heterogeneity that is concomitant to the glass transition. In the past, dynamic heterogeneity of glass forming liquids has been investigated by many authors [13][14][15][16][17][18][19][20][21][22] using two point correlation functions such as, intermediate scattering function, mean-squared displacements, non-Gaussian parameters etc. Many investigations have argued that conventional two point correlation functions are not adequate for fully describing the spatially heterogeneous dynamics, and the relation between length and time scales of relaxation in supercooled and glassy liquids [23]. To measure the extent of dynamic heterogeneity, several authors have employed four point correlation functions in computational studies of model glass formers [14,17,[24][25][26][27][28][29][30][31][32][33][34][35].
There are large variety of materials show glass transition, including polymers [36]. This work focus on the glass transition in polymer melts, which are well known glass formers and have numerous commercial, industrial, and daily life applications [37]. As polymers have many degrees of freedom, it is computationally expensive to address glass transition in polymer systems. One of the simplest model of a polymer liquid is a system of linear chains. We have performed equilibrium molecular dynamics(MD) simulations of supercooled system of linear polymer melt of a model polymer proposed by Grest and Kremer(GK) [38], which was used for earlier studies of glass transition by Bennemann et al. [39]. They studied the divergence of α-relaxation time at temperatures above T c . Increasingly correlated motion, which was more pronounced well above T c , has been studied by Bennemann et al. [40] in the same model of polymer melts but at higher densities. Spatially correlated dynamics and string like motion of the same system has been studied by Glotzer and coworkers [21,41]. Starr et al. [42] compared string lengths to the length scale of Adam-Gibbs theory [43] and random first order transition theory [44][45][46] of the glass transition. Extensive MD simulations has been performed by Leporini et al. [47][48][49][50][51] to measure the spatial correlated motion of the particles in the supercooled polymer melt. Jump like motions near the T c are investigated by Helfferich et al. [52,53] in a glassy polymer melt using continuous time random walk approach. It is now interesting to extend earlier studies to further establish the connection between microscopic structure and dynamics of this system especially by calculating the four-point correlation functions. The concept of four point correlation functions has been proposed by Dasgupta et al. [54] in their simulation studies. Donati et al. [55] showed that the growth of correlation length can be given by the volume integral of four point correlation functions. Four point density correlation functions have been computed to understand increasingly correlated motion and dynamic correlation length in various glass forming liquids [14,34,35,42].
We describe basic theory of correlation functions in Sec.2; the model and simulations detail is presented in Sec.3; finally we will present the results and analysis part of this work in Sec.4.

Basic theory of static and dynamic correlations
For homogeneous and isotropic systems, radial distribution function g(r) is only function of separation of a pair of particles, i.e., r and then the g(r) is defined as [56] The g(r) of polymer melt system can be expressed as g(r) = g intra (r) + g inter (r), where sum for g intra (r) takes contribution from monomers of same chain, and sum for g inter (r) takes contribution from monomers that belong to different chains [57]. The single particle relaxation dynamics of the system can be analyzed using incoherent intermediate scattering function [56], which can be defined as and the structural α-relaxation time, i.e., τ α can be defined as τ α = ∞ 0 F s (k, t) dt. The structural dynamic susceptibility χ 4 s (k, t) can be defined from the variance of F s (k, t) [27,31,58]  We can extract the length and time scale of dynamic heterogeneity from χ 4 s (k, t) of a model glass-former.

Model and simulation details
We simulated coarse-grained bead-spring polymer melt model system for modeling the polymer chains, each linear chain consists of 10 monomers. The non-bonding interactions are modeled by Lennard-Jones(LJ) potential [59] that is given by The LJ potential is truncated at distance of r cut = 2 × 2 1/6 σ, where σ and are diameter of a particle and depth of the potential respectively. The bonding interactions are modeled by finitely extensible non-linear elastic(FENE) potential [38] which is given by , where k 0 = 30 σ −2 and R = 1.5σ. We integrated the Newton's equations of motion using Velocity-Verlet algorithm [59] with time step δt = 0.002. We simulated a linear polymer melt system of N c = 1000 polymer chains. The equilibration time used in the simulation is the time required for relaxation of end-to-end vector to 10% of its initial value at all temperatures [52]. The system is studied at the temperatures T = 2.0, 1.0, 0.90, 0.80, 0.70, 0.60, 0.50, 0.45, and 0.40 at fixed density ρ = 0.85 in NVE ensemble. All the parameters reported from our simulation study are in LJ simulation units, i.e., the physical quantities are represented as follows: distance in units of σ; the reduced density, ρ * = ρσ 3 ; temperature, Fig.1(a) shows the g(r) of the system at various temperatures, which shows a broad peak at r = 1.0 at high temperatures. When temperature is reduced, the first peak shows prominent splitting from temperature T = 1.0. The first peak at low temperatures in Fig.1(a) corresponds to average bond length(r bond ≈ 0.96) and second peak is due to LJ interaction this is around r = 1. When temperature decreases towards T c , first peak and the subsequent peaks grow that shows formation of the secondary structure in this system. As polymers can form complex structures within same molecule, it is meaningful to look separately in structures formed within each polymer molecules and also the structure formed by surrounding molecules. In the inset of Fig.1(b), g intra (r) shows the structure formed around a bead of a chain by monomers of the same chain. Like g(r), the peaks in g intra (r) grow while reducing the temperature towards T c . The average radius of gyration of the polymer melt chains is R g ≈ 1.44 in equilibrium, which has weak temperature dependence. The second peak and subsequent peaks grow at low temperatures due to the confinement of a chain via LJ interaction among monomers of the same chain, that results in the globule like structure of chains in these simulations. The peak of the g intra (r) near r = 2 also broadens showing multiple secondary structures within a chain similar to binary LJ glass forming systems [60] near T c . In Fig.1(b), we show g inter (r), which shows the structure formed by the monomers belong to neighboring chain that surrounds a tagged chain. It is clear from this plot that g inter (r) contributes more to the g(r) in comparison with g intra (r) in the formation of secondary structure. In this case also the peaks that is around r = 2 broadens to show formation of the secondary structure. Now it is interesting to look into how the relaxation dynamics affected by these multiple secondary structures that is formed at low temperatures.  Fig.1(b) shows g inter (r) and inset of Fig.1(b) shows g intra (r).

Results and discussions
We have analyzed single particle relaxations by computing the F s (k, t), from Eq.2, which is plotted in Fig.2. The wavenumber used in the calculation is k = 6.99, this corresponds to the maximum of static structure factor S(k). The F s (k, t) fits well with exponential function at higher temperature T = 4.0, below this temperature, it is fitted by empirical Kohlrausch-Williams-Watts function of the form f (t) KW W ∝ exp(−( t τ ) β ) with value of exponent β varying from 0.87 to 0.61, as temperature decreases for this polymer melt system. At temperatures, above and near T c , particles shows pronounced two step relaxations, showing β and α regime of the relaxation function in supercooled linear polymer melt system. The onset of plateau in the F s (k, t) starts at T = 0.50, which signifies the heterogeneous relaxation dynamics in the system. We computed τ α at aforementioned temperatures; the τ α vs T curve is fitted using MCT and Vogel-Fulcher-Tammann relations, i.e., τ α ∝ (T − T c ) −γ and τ α ∝ exp(DT 0 /(T − T 0 )). The fitted parameters for the linear polymer melt system, at given thermodynamical conditions, of these two equations are given as follows: MCT temperature, T c ≈ 0.40; the MCT exponent, γ ≈ 1.65; the fragility index, D ≈ 2.34; the dynamic divergence temperature, T 0 ≈ 0.30. Many parameters we obtained are comparable to the values obtained in the work of Bennemann et al. [39] in NVT and NPT ensembles, however, our calculations are for a larger system and in NVE ensemble. The ideal glass transition temperature T g and extrapolated dynamic divergence temperature T 0 are related through the relation T g ≈ 1.2T 0 [42]. The approximate T g for this system is 0.36. Now we look into the four-point correlation functions to understand the increasingly correlated motion of particles in this linear polymer melt. In Fig.3, we have given self part of four point χ 4 s (t) of monomers of linear polymer chains. The height of the peak and corresponding time grows in χ 4 s (t), while reducing temperature towards the T c in this system. Interestingly, time corresponding to the peak of χ 4 s (k, t) in our system shows same variation as the α-relaxation time. Berthier et al. [61] reported experimentally that the time corresponding to the peak of χ 4 (t) is equivalent to the average relaxation time of the system, i.e., t 4 max ≈ τ α , where t 4 max time corresponding to the maximum in χ 4 s (t), in supercooled glycerol and colloidal hard spheres. The pronounced peak of χ 4 s (k, t) reveals the increasingly correlated motion of the particles in supercooled polymer melt, which is due to the spatially heterogeneous dynamics. This increasing height of the peak of χ 4 s (k, t) towards T c gives the signature of the growing correlation length in the system. This growing correlation length and time quantify the spatio-temporal heterogeneity of the system in the vicinity of T c .

Conclusion and summary
We have reported results of the long time molecular simulation study of linear polymer melt system of 10 3 particles in micro-canonical ensemble near its glass transition temperature. Our studies looked into the microscopic structure of the slowly relaxing domains near glass transition temperature T c . Explicit analysis of the radial distribution functions of linear polymer melts shows formation of multiple structures in this system as glass transition is approached. We have computed the structures formed by intra and inter chain monomers using g intra (r) and g inter (r). Analysis of g intra (r) and g inter (r) revealed the role played by these monomers in formation of secondary structures. As reported in earlier studies [39], we also find appearance of plateau region in the F s (k, t) and divergence of α-relaxation time [39] in linear polymer melts at low temperatures, which is a signature of spatially heterogeneous dynamics in the system. To further characterize spatio-temporal dynamic heterogeneity, in this system of polymer melt, we have computed structural dynamic susceptibility χ s 4 (k, t). The calculated χ s 4 (k, t) gives the correlation volume of a system, which diverges at low temperatures, signifies the increasingly correlated motion of the particles in this linear polymer melt system.

Acknowledgments
We acknowledge central HPC facility at IIT Mandi for computational support.