Dynamics and diffusive–conformational coupling in polymer bulk samples and surfaces: a molecular dynamics study

The impact of free surfaces on the mobility and conformational fluctuations of model polymer chains is investigated with the help of classical molecular dynamics simulations over a broad temperature range. Below a critical temperature, T*, similar to the critical temperature of the mode coupling theory, the center-of-mass displacements and temporal fluctuations of the radius of gyration of individual chains—as a fingerprint of structural reconfigurations—reveal a strong enhancement close to surfaces, while this effect diminishes with increasing temperature and observation time. Interpreting conformational fluctuations as a random walk in conformational space, identical activation enthalpies for structural reconfigurations and diffusion are obtained within the error bars in the bulk and at the surfaces, thus indicating a coupling of diffusive and conformational dynamics.


Simulation techniques
Classical MD simulations were performed 5 on a model polymer glass using a bead-spring model of the Kremer-Grest type [4]. A truncated and shifted Lennard-Jones (LJ) potential acts between all particles, where the constant c guarantees continuity for all r , and r c denotes the truncation distance with r c = 2 r min = 2 × 2 1/6 σ . Along a chain, adjacent beads further interact through a finite extensible nonlinear elastic (FENE) potential, ensuring the connectivity of a chain: with k = 30 /σ 2 and R 0 = 1.5σ . These parameters avoid chain crossing and high-frequency modes [9], thus allowing for simulations of a large system over sufficiently long periods of time. In the following, all results are given in reduced LJ units, where energy, length, mass and k B are set to unity, leading to a unit of time of t * = σ (m/ ) 1/2 . All simulations were performed in the isobaric-isothermal ensemble (NpT) using a Nosé-Hoover thermostat and barostat [27,28] with a time step dt = 0.005. Starting configurations were obtained using a method proposed by Kremer and Grest [4]; 1730 chains of 30 beads, each with a fixed bond length, were randomly placed in a cubic box with periodic boundary conditions applied in the xand y-direction and two surfaces to the vacuum in the z-direction in order to simulate a free-standing thin film. Overlaps were removed by a short MD simulation run at pressure P = 1 with a softer potential U (r ) = A [1 + cos(π r/2 1/6 σ )] (2 1/6 σ − r ), where (r ) denotes the Heaviside function. After that pressure control was disabled and the system was equilibrated at a high temperature T = 1 for t = 6 × 10 6 dt, cooled to T = 0.4 within t = 1 × 10 5 dt, where the system further relaxed for another t = 6 × 10 6 dt, before it was cooled again to T = 0.2 in t = 1 × 10 5 dt to approach a low-temperature starting point. During cooling the volume of the cell clearly shows the signatures of a glass transition. To determine T g , the system was cooled from T = 1 to 0.1 with different cooling rates between 0.9/10 4 dt and 0.9/10 6 dt, leading to 0.42 T g 0.46. The absence of long-range order, which can be visualized by evaluating the pair distribution function, additionally corroborates amorphicity. After the first equilibration the system was heated from T = 0.2 to the desired temperature, where the system was relaxed again for t = 8 × 10 6 dt, i.e. until the potential energy ceased to show further reduction as a function of time. Afterwards production runs were performed for t = 1 × 10 8 dt, making snapshots of the configurations every 1 × 10 6 dt, while the total system enthalpy remained constant within the magnitude of typical thermal fluctuations. While the present study encompasses 17 temperatures in the interval 0.39 T 0.62, further tests were performed to exclude finite size effects by employing twice as many chains, leading to similar results. Effective diffusion constants D cm (z) as a function of the initial centerof-mass distance from the film center, z, at different sample temperatures. At lower temperatures, a pronounced increase appears near the surface (z ≈ ±15), which decreases as temperature is increased. Nevertheless, mobility generally increases with rising temperature. The solid lines are guides to the eye.

Polymer diffusion
The mobility of a chain i can be characterized by the mean-squared displacement (MSD) of its center-of-mass (cm-given by R cm,i ), x 2 cm,i (t) = ( R cm,i (t) − R cm,i (t 0 )) 2 , while we set the starting point, t 0 , of our production runs (which directly follow sufficient relaxation, as described above) equal to zero for convenience. In principle, one could also consider other MSDs, e.g. the MSD of monomers; however, the latter reach the diffusive limit later than the centerof-mass MSD due to an evolving intermediate subdiffusive regime [13]. As a quantitative measure of mobility, an effective long-time limit diffusion constant, D i , of chain i can be calculated throughout the production runs, i.e. D i = x 2 cm (t)/6t for t 0. Averaging D i with respect to the initial distance of the chain's center-of-mass from the film center, z cm,i , we obtain which is shown in figure 1 for representative temperatures ranging from T = 0.39 (well below T g ) up to T = 0.62. At low temperatures, we observe a strong increase of effective diffusion near the surfaces, which diminishes-in penetration depth as well as magnitude-at higher temperatures, while a small surface enhancement seems to remain. To quantify this point, figure 2 shows an Arrhenius plot of the effective diffusion constant in the core region of the cell, D cm,b = D cm,|z| 5 (which, in fact, agrees well with the bulk diffusion constant determined from simulations with periodic boundary conditions in all directions), and two surface regions 6 , D cm,s = D cm,|z| (max(z)−2) , respectively. Toward lower temperatures, bulk and surface values start to deviate notably; the ratio of D cm,s and D cm,b decreases from D cm,s /D cm,b ≈ 45 at T = 0.39 to D cm,s /D cm,b ≈ 1.3 at T 0.5. The complete temperature dependence is depicted in an Arrhenius plot (figure 2), which clearly reveals-for bulk and surface diffusion-a transition from the high-temperature liquid regime to glassy dynamics at a critical temperature, T * D ≈ 0.47. Larger than T g , the latter has great similarities with the critical temperature of the mode coupling theory [29], T C (see e.g. [30] for the experimental metallic glass equivalent). From the slopes, the corresponding enthalpies of activation, H D , are readily calculated, as included in figure 2. In the glassy state, H D is generally much higher in the bulk ( H D ≈ 5.0(3)) than in the surface region ( H D ≈ 2.3(2)). In the high-temperature regime, this difference diminishes to H D ≈ 1.6(1) and H D ≈ 2.0(1), respectively. Generally, the increased mobility in surface proximity is in good agreement with recent experiments [24,25] and computer simulations [18,19]. The dependence of the diffusion profile on surface distance has some similarities to the findings in metallic glasses [20,21,31], which also showed a remarkable increase when compared with the bulk, thus indicating the basic universality of this scenario. Enhanced surface mobility decreases and seems to approach a lower threshold above T ≈ 0.5. This can be attributed to the intermixing at high temperatures, whereby chains do not stay near the surface but can move into the cell and vice versa. It should be mentioned that in our studies the open surface has only minor effects on the corresponding density, which raises questions about the applicability of free volume concepts [32] as an explanation. Even though the behaviors of bulk and surfaces are clearly different, there is a continuous transition between bulk and surface properties, while indications for different phases in the bulk and at the surface are absent.

Structural reconfigurations
The radius of gyration, R G (t), is a common measure to describe the size and conformation, viz. dimension, of a chain at a particular time, t. It can readily be calculated from the monomer with time always belong to conformational changes and thus structural reconfiguration. Fluctuations of R 2 G are thus expected to reach a maximum, when time and/or temperature are large enough to allow for a most complete sampling of conformational phase space. The liquid state above the critical temperature of the mode coupling theory, T C , by definition always allows for such a complete sampling, while lower temperatures successively require longer observation times, until full sampling ceases on the observation timescales in the glassy state (i.e. for T T g ). We calculate the expectation value of fluctuations of the radius of gyration of chain i at time t by evaluating Here, R 2 G,i t denotes the time-averaged radius of gyration of chain i, while a number of n snapshots are taken into account within the observation time t. Applying the same procedure with respect to the initial z-position yields the fluctuation depth profile, which is shown indicates the equivalence of observation time and temperature-staying longer at a specific temperature has the same effect as shorter observation times at higher temperatures [34], and vice versa. It is worth noting, though, that differences in the fluctuation amplitude between bulk and surface are apparent at low temperatures and short observation times only, while transitions with increasing temperature and observation time occur smoothly. It should be mentioned that at sufficiently long observation times, a chain samples the whole system and thus all possible distances from the surface, which is probably the reason for the merging of bulk and surface properties. The mean radius of gyration of the chains, which have a length of 30 in this study, is R 2 G ≈ 7. This value is only slightly dependent on temperature.
A corresponding analysis has been performed with the end-to-end distance of the chains (not shown here) leading to qualitatively similar results. Figure 4 shows the temporal evolution of the average fluctuations in the bulk and at the surfaces, respectively, for different temperatures. At all temperatures, fluctuations increase with time to a maximum level of ≈3, while the timescale of convergence strongly increases with reduced temperature, until-at the lowest temperatures-saturation cannot be reached within the simulation time. For the sake of clarity, only representative temperatures are shown here; at higher temperatures, the values conglomerate near the T = 0.62 ordinate. Below T ≈ 0.48 surface fluctuations are always higher than in the bulk, whereas at higher temperatures these differences level.
We proceed by a more detailed analysis of the shape of the curves of R 2 G (t) in figure 4. As indicated by the solid lines, the data points are well fitted by a stretched exponential Kohlrausch-type [33] function of the shape where the stretching exponent β is determined to be ≈0.44 (5). A very similar stretchedexponential behavior has been reported recently for the temporal evolution of the tensor of elasticity in proximity to the glass transition temperature [35], and can be rationalized based on the kinetics of phase space sampling. While we will come back to this aspect later, we first focus on the relaxation times, τ , determined by fitting equation (5) to the data in figure 4. Cast in an Arrhenius plot (figure 5), the latter reveal a temperature behavior very much akin to the diffusion constants in figure 2. Clearly, the temperature dependence can be divided into two regimes-above and below a critical temperature, T * R , respectively-while in each of the regimes an Arrhenius behavior with characteristic activation enthalpies, H R , for the bulk solid and the surfaces, respectively, prevails. Qualitatively, the activation enthalpies for conformational changes on the surfaces prove to be significantly lower than in the bulk, which is physically reasonable due to reduced constraints imposed on surface polymers, when compared with the bulk. Also, the occurrence of a critical temperature, T * R , at which polymer kinetics significantly change is qualitatively understandable, when relating T * R to the critical temperature of the mode coupling theory, T C . Above T C the matrix surrounding of an individual polymer appears as liquid on the timescale of MD, and thus imposes only a reduced constraint on to a conformational change (as reflected by a reduced activation enthalpy H R ), when compared to below T C .

Relation between diffusion and structural relaxation
A striking feature of our diffusion and structural reconfiguration data is the fact that the Arrhenius plots of the corresponding relaxation times 7 (figures 2 and 5) reveal identical activation enthalpies-above T C , below T C , in the bulk and at the surface. This constitutes a strong indication that the same molecular relaxation process underlies both scenarios from the glassy up to the liquid state. In fact, coupling between different relaxation processes in glasses has recently been discussed based on the interpretation of experimental data [36,37], 9 and is-after all-one of the foundations of the mode coupling theory [29]. A similar coupling of diffusion constant and the orientational relaxation time of end-to-end vectors (which is on a similar length scale to the radius of gyration) in the liquid regime above T C has been reported [38], confirming the close connection between conformational dynamics and diffusion. Nevertheless, it has also been reported that the decay time of the incoherent scattering function, which is interpreted as a structural relaxation time on the length scale of a monomer, has a different temperature dependence. Based on our present results, we interpret the coupling of diffusion and conformational relaxation as follows. Individual polymers basically sample their conformations (and thus conformationally relax) by a random walk in conformational space (as corroborated by β ≈ 0.5, as in [35]), while each conformational change is associated with a random walk in real space. The latter is, in fact, basically one of the foundations of the well-established reptation model [6].

Conclusion
To summarize, MD simulations on the dynamics of bulk samples and thin films of a model polymer system were investigated with regard to diffusion and conformational relaxation. Enhanced diffusivity and structural reconfigurability were detected in surface proximity, when compared with the bulk. Evaluation within Arrhenius plots revealed a transition in both dynamic properties at a temperature identified with the critical temperature of the mode coupling theory, T C . Identical activation enthalpies in all temperature regimes within the error bars were interpreted as indicating a common underlying dynamics. The latter was rationalized as a random walk in conformational space.