The glass transition temperature of confined and free-standing polymer films of varying thickness is studied by extended molecular dynamics simulations of bead–spring chains. The results are connected to the statistical properties of the polymers in the films, where the chain lengths range from short, unentangled to highly entangled. For confined films, perfect scaling of the thickness-dependent end-to-end distance and radius of gyrations normalized to their bulk values in the directions parallel and perpendicular to the surfaces is obtained. In particular, the reduced end-to-end distance in the perpendicular direction is very well described by an extended Silberberg model. For bulk polymer melts, the relation between the chain length and Tg follows the Fox–Flory equation. For films, no further confinement induced chain length effect is observed. Tg decreases and is well described by Keddie’s formula, where the reduction is more pronounced for free-standing films. It is shown that Tg begins to deviate from bulk Tg at the characteristic film thickness, where the average bond orientation becomes anisotropic and the entanglement density decreases.
In a majority of applications (commodity), amorphous polymeric materials are in the glassy state. Because of that, the glass transition temperature Tg region is of crucial importance,1 and thus, fundamental and applied research on the glass transition, in general, and of polymers, in particular,2–4 remains an important field of study in condensed matter physics. Furthermore, when an amorphous polymer in the liquid state is cooled toward Tg, the viscosity dramatically increases in a non-Arrhenius way.5–7 The variable polymer chain length provides an additional parameter for systematic studies compared to low molecular weight glass formers.8,9 The special role polymers can play for a detailed view on the glass transition is illustrated in a recent comprehensive study of the temperature-dependent relaxation dynamics of different polymers on chain length, chain flexibility, and chemical properties.9 Similarly, Xu et al.10 employed extensive molecular dynamics (MD) simulations based on the model of Ref. 11 to investigate the role of chain stiffness on the glass transition. Experimentally, typical ways to determine Tg of polymers are by differential scanning calorimetry (DSC)12 or by thermomechanical analysis (TMA),13 which can lead to deviating results.9 This holds for both simple, i.e., low molecular weight, glass formers and polymers. However, despite this huge research effort, the nature of the glass transition is still not fully understood.4,14–22 These studies mainly refer to bulk systems. Free-standing or supported thin films pose even more challenges and have led to—partially—contradictory claims in the literature.23–27 While mostly it is reported that Tg is reduced for very thin films, in some cases, depending on chemistry and substrate, even an increase was observed.25,28–30 It is generally understood and stressed by the previous studies that this is due to the relevance of sample preparation and surface coupling or even environmental conditions. This leads to questions concerning comparability and equilibration of such submicron thick films. Due to confinement and surface interaction, a significant impact on dynamic and structural polymer properties is observed. For instance, in-plain chain extensions are only weakly affected by confinement, while in the perpendicular direction, the chain extension gradually reduces from the bulk value with a decreasing film thickness.17,31–36 This also leads to chain mobility modifications,35,36 which depend on the distance to the surface. Especially for very long, highly entangled polymers, these problems become very complex. The aim of this Communication is to systematically investigate the dependence of Tg of confined and free-standing films on chain length and film thickness.
Taking all these interconnected problems, computer simulations of thin polymer films offer insight under perfectly controlled conditions. Full time dependent coordinates of all particles are available, allowing us to investigate the structure and molecular motion (viscosity) and the glass transition under a variety of situations.10,11,16,37,38 So far, however, most computational studies on polymer films have focused on short and unentangled chains38 because the computing time for equilibration increases dramatically with chain length and systems’ complexity increases.
We here apply a recently developed efficient hierarchical methodology to equilibrate the highly entangled melts of long polymer chains in bulk,39,40 and confined and free-standing polymer films.41 The required computer time scales linearly with system size, independent of chain length. While not restricted to standard bead–spring chains, we here employ this approach to weakly semiflexible bead–spring chains,42,43 where the entanglement length Ne = 28 beads at the standard melt density of ρ0 = 0.85σ−3. The reference temperature for equilibration and cooling is T = 1ɛ/kB. Throughout the whole paper, the Lennard-Jones (LJ) units of length (σ), energy (ɛ), and time (τ) are used.
For equilibration and structural properties at T = 1ɛ/kB, we use the mentioned standard semiflexible bead–spring model. This model, here referred to as model I, is widely used in the literature, and its properties are well documented.17,42,44–51 Unfortunately, this cannot be used to study free surfaces and the glass transition. Model I only contains repulsive non-bonded interactions and displays an atypical chain stretching upon cooling. To correct for that, we employ a variant, model II,11 with attractive non-bonded interactions, allowing for free surfaces where the pressure is kept at52–57, P = 0.0ɛ/σ3 at least along the direction perpendicular to the interfaces of films, and a modified bond angular potential, ensuring that the chain conformations only very weakly depend on T. The details of the parameterization are chosen such that at T = 1ɛ/kB and ρ = 0.85σ−3, the chain conformations and bead packing between the two models are indistinguishable. This close resemblance to the standard semiflexible bead–spring model allows us to switch “on the fly” between model I and model II. We can, thus, take advantage of the huge body of available data of model I. A detailed comparison of the models is given in the supplementary material. Model II captures the major features of glass-forming polymers, e.g., that the viscosity and relaxation time dramatically increase in a non-Arrhenius way close to Tg11,58 as it is characteristic for fragile polymeric glass formers. For confined polymer films, two confining planar, structureless repulsive walls with a 10-4 LJ potential41,59,60 are introduced. This is sufficient to investigate the generic conformational properties of the films.
The ESPResSo++ package61,62 is used to perform MD simulations with a Langevin thermostat in the NVT and NPT ensemble, with a Hoover barostat for the latter. Our polymer melts contain nc chains of N monomers (ncN = 105 for N ≤ 100 and nc = 1000 for N = 500 and N = 2000), ranging from unentangled (N < Ne) to highly entangled (N ≫ Ne) films. The film thicknesses h range from thick to very thin films, being the root mean square (rms) radius of gyration for bulk melt chains. The smallest h ≈ 9.0σ < 2dT (see Fig. 1), dT = 5.02σ being the tube diameter.43 In all cases, h is measured along the z direction according to the concept of the Gibbs dividing surface that has been applied to identify the interface between two different phases,41,63 while the periodic boundary conditions are applied in the x and y directions. All confined polymer melts of short chains N ≤ 100 could easily be generated via a brute-force equilibration as in the bulk.42
Snapshot of a configuration of fully equilibrated free-standing film (nc = 1000, N = 2000 ≈ 72Ne) of at T = 1.0ɛ/kB. The dashed box is shown for better visualization.
Snapshot of a configuration of fully equilibrated free-standing film (nc = 1000, N = 2000 ≈ 72Ne) of at T = 1.0ɛ/kB. The dashed box is shown for better visualization.
We first analyze the chain conformations for confined melts at T = 1.0ɛ/kB in comparison with the bulk properties as shown in Fig. 2. For each component of the rms end-to-end distance and radius of gyration, we find excellent data collapse onto a universal master curve. This is an additional indication that all confined polymer films are, indeed, fully equilibrated. While Re,‖(N, h) and Rg,‖(N, h) remain unchanged almost down to , this is different for the perpendicular component. There, we observe a gradual decrease already for much thicker films originating from chains close to the surface. Around , this turns into a linear decrease with h for Re,⊥(N, h) (Rg,⊥(N, h)), while weakly increases for Re,‖(N, h) (Rg,‖(N, h)).33 Earlier, in 1982, Silberberg31 argued that chains next to a wall can be treated as unperturbed random walks located with their center at the wall and folded back on one side of the wall. This changes Re,⊥(N, h) described by the scaling function ; see the supplementary material; however, it does not affect the parallel component. This idealized picture works well in the thick film regime down to . Sussman et al.35 extended Silberberg’s hypothesis from one wall to two walls, keeping the reflecting boundary conditions in the perpendicular direction and the Gaussian weight factor for each contribution. Their extended theoretical prediction of , shown in the supplementary material up to the second-order correction term, is supported by our data covering the thick and thin film regimes. For ideal chains, , while the distribution of is not exactly Gaussian anymore. Therefore, our data on Rg,⊥(N, h) slightly deviate from . Furthermore, a detailed analysis of our data in the parallel direction reveals that the scaling behavior of Re,‖(N, h) and Rg,‖(N, h) is better described by f‖(x) = 1 + chx2 instead of f‖(x) = 1 + chx, with as in Ref. 64. Here, (see Fig. S1 in the supplementary material) is the excluded volume screening length. is the degree of interdigitation of different polymers, also called generalized polymerization index64 in bulk melts.
(a) and (b) Two components of rms end-to-end distance and radius of gyration rescaled to the bulk value, (a) and (b) in the directions parallel (α = ‖) and perpendicular (α = ⊥) to the walls, plotted vs h/ξ(0)(N) and , respectively. (c) Orientational order parameter S and reduced number of kinks, , plotted vs h. In (a) and (b), the theoretical predictions f‖(x = h/ξ(0)(N)) = 1 + chx2 with the fitting parameter ch ≈ 0.22, ,31 and 35 (cf. text) are shown by the curves for comparison. In (c), , 22.28(1.31), and 91.78(2.92) for N = 100, 500, and 2000, respectively, and the dashed lines are drawn to guide the eye. All data are for T = 1.0ɛ/kB.
(a) and (b) Two components of rms end-to-end distance and radius of gyration rescaled to the bulk value, (a) and (b) in the directions parallel (α = ‖) and perpendicular (α = ⊥) to the walls, plotted vs h/ξ(0)(N) and , respectively. (c) Orientational order parameter S and reduced number of kinks, , plotted vs h. In (a) and (b), the theoretical predictions f‖(x = h/ξ(0)(N)) = 1 + chx2 with the fitting parameter ch ≈ 0.22, ,31 and 35 (cf. text) are shown by the curves for comparison. In (c), , 22.28(1.31), and 91.78(2.92) for N = 100, 500, and 2000, respectively, and the dashed lines are drawn to guide the eye. All data are for T = 1.0ɛ/kB.
This change of conformation of confined films also affects the orientational bond distribution in the chains as analyzed by the bond orientational order parameter S = (3⟨cos ϕ⟩ − 1)/2, where ϕ is the angle between each bond and the z axis. Furthermore, entanglements will be affected by the change in chain self-density. Applying a primitive path analysis65 to obtain entanglement points (significant kinks)66,67 along the primitive paths (PPs) reveals this. The results of S and the average number of significant kinks Nkink as a function of h are shown in Fig. 2(c). S is independent of N for a given h, indicating that the local packing is not affected by N. In contrast, the estimates of Nkink for N = 500 and 2000 follow a master curve and deviate strongly from that of shorter chains of only a few entanglement lengths (N = 100 ≈ 3.6Ne). However, both datasets deviate from the bulk below and around h = hc ≈ 20.0σ. Obviously, hc is not related to , differently from the prediction in Ref. 35. The chains are less strongly oriented on the scale of bond vectors between monomers compared to the scale of the end-to-end vector.36
To model free-standing films in vacuum and to study the glass transition of polymer films, we apply model II to all confined polymer films shown in Fig. 2 at T = 1.0ɛ/kB. After a short initialization in the NVT ensemble, the confined polymer films are further relaxed for about 13τe in an NPT ensemble with a fixed wall distance at pressure P = 0.0ɛ/σ3. This led to a marginal adjustment of the lateral extensions of the films. Then, free-standing films are obtained just by removing the wall potential while keeping their lateral dimensions fixed (for more details, see Ref. 41). After that, we perform MD simulations of confined films and bulk melts in the NPT ensemble and those of free-standing films in the NVT ensemble, keeping the dimensions of films constant. The component of pressure tensor along the direction perpendicular to the interfaces Pzz for all films fluctuates around zero.
To study the glass transition, we follow exactly the very same cooling protocol as for our previous bulk studies of the same polymer model.11,68,69 We apply stepwise cooling,54 which results in a cooling rate Γ = ΔT/Δt = 8.3 × 10−7ɛ/(kBτ). The temperature is reduced in steps of ΔT = 0.025ɛ/kB from T = 1.0ɛ/kB to 0.2ɛ/kB with a relaxation time between each step of Δt = 30 000τ ≈ 13τe ( being the entanglement time defined by the tube model70 and reptation theory,71,72 and τe ≈ 2266τ estimated from simulations42,43), i.e., subchains of the order of Ne can relax easily at higher temperatures close to T = 1.0ɛ/kB.
To determine Tg, we perform a hyperbolic fit74 to the density ρ(T) change with temperature , where c, T0, a, b, and f are the fitting parameters. The density is defined by ρ(T) = ncN/(h(T)L(T))2, where L(T) is the lateral dimensions and h(T) is the effective film thickness at the temperature T; see the supplementary material. Adopting this fit, Tg is defined by either Tg = T0 or the intersection point of two tangents drawn at the high and low temperatures. Both give the same estimate within fluctuations for all systems studied here; see, for example, Figs. 4(a) and 4(b) for confined and free-standing films of two selected film thicknesses h for the longest polymer chains of N = 2000. For comparison, the estimates of of the corresponding bulk melts11,68,69 are shown in Fig. 3. As in experiment and other simulations, our data are well described by the Fox–Flory relation.9,38,73 decreases with N for N < 2Ne, while for N ≳ 2Ne.
Estimates of from the density ρ(T) change, plotted as a function of N at P = 0ɛ/σ3 for polymer melts in bulk. The Fox–Flory relation73 with K = 0.579(59) and is shown by a black curve for comparison. The uncertainty of is indicated by a gray shaded region.
Estimates of from the density ρ(T) change, plotted as a function of N at P = 0ɛ/σ3 for polymer melts in bulk. The Fox–Flory relation73 with K = 0.579(59) and is shown by a black curve for comparison. The uncertainty of is indicated by a gray shaded region.
Figures 4(c) and 4(d) show the results of Tg(N, h) of confined and free-standing films depending on h and N following the same data analysis as shown in Figs. 4(a) and 4(b). Obviously, Tg(N, h) covering the range from unentangled to highly entangled chains only very weakly depend on N, while the reduction with decreasing h in both confined and free-standing films is clearly observed. The data are well described by30, with an exponent δ = 2.0(1), consistent with δ = 1.8(2) for the supported PS films.30,75 The characteristic length h0 for the free-standing films is about 1.5 times that for the confined films due to the different chain mobility near the surface.19, h0 also depends on the chemical composition of polymer chains.76 The deviation of Tg(N, h) from around h ≈ 20σ fits well to the observed changes of the bond orientation and the reduction in entanglements in Fig. 2(c). Most notably, the relative depletion of Tg(N, h) seems to be almost independent of chains’ length as observed from experimental studies,76 unlike the value of Tg(N) in bulk. Alternatively, one can also determine Tg from the bilinear fits of the total potential energy Utot(T) change54,77 (see the supplementary material). The estimates of Tg for N = 2000 are shown in Fig. 4. The Tg reduction remains the same; however, the uncertainty is much larger.
(a) and (b) Rescaled monomer density ρ(T)/ρ(T = 1) in confined (a) and free-standing polymer films (b) of h/R(0) (N) ≈ 4.3 and 0.3, plotted vs T for N = 2000. The hyperbolic and tangent fits are represented by the curves and lines, respectively. The estimates of Tg(N, h) via the hyperbolic fit are indicated by the arrows including the uncertainty. In the insets of (a) and (b), full data for 0.2 ≤ kBT/ɛ ≤ 1.0 are shown. (c) and (d) Tg(N, h) rescaled by (see Fig. 3) including the error bars, plotted as a function of h for confined (c) and free-standing (d) films. The uncertainty of is indicated by a gray shaded region. The formula proposed by Keddie et al.30 with δ = 2.0(1) and h0 = 1.67(23)σ in (c) and h0 = 2.55(35)σ in (d) are shown by the dashed and solid curves, respectively. The estimates of Tg(N, h) from the potential energy Utot(T) change for N = 2000 are also shown in (c) and (d) for comparison. At T = 1.0ɛ/kB, σ3ρ(T) ≈ 0.854 and 0.858 for h/σ ≈ 130.0 and 9.0, respectively, in (a), and σ3ρ(T) ≈ 0.853 and 0.852 for h/σ ≈ 130.3 and 9.0, respectively, in (b).
(a) and (b) Rescaled monomer density ρ(T)/ρ(T = 1) in confined (a) and free-standing polymer films (b) of h/R(0) (N) ≈ 4.3 and 0.3, plotted vs T for N = 2000. The hyperbolic and tangent fits are represented by the curves and lines, respectively. The estimates of Tg(N, h) via the hyperbolic fit are indicated by the arrows including the uncertainty. In the insets of (a) and (b), full data for 0.2 ≤ kBT/ɛ ≤ 1.0 are shown. (c) and (d) Tg(N, h) rescaled by (see Fig. 3) including the error bars, plotted as a function of h for confined (c) and free-standing (d) films. The uncertainty of is indicated by a gray shaded region. The formula proposed by Keddie et al.30 with δ = 2.0(1) and h0 = 1.67(23)σ in (c) and h0 = 2.55(35)σ in (d) are shown by the dashed and solid curves, respectively. The estimates of Tg(N, h) from the potential energy Utot(T) change for N = 2000 are also shown in (c) and (d) for comparison. At T = 1.0ɛ/kB, σ3ρ(T) ≈ 0.854 and 0.858 for h/σ ≈ 130.0 and 9.0, respectively, in (a), and σ3ρ(T) ≈ 0.853 and 0.852 for h/σ ≈ 130.3 and 9.0, respectively, in (b).
In summary, we confirm that the dependence of Tg of bulk polymer melts on the chain length N in the range from unentangled to highly entangled can be well described by the Fox–Flory equation. The scaling prediction of the chain extensions in the directions parallel and perpendicular to the surface of confined films is verified in both thick and thin film regimes. We show that Tg, S, and Nkink (for films of highly entangled polymer melts) start to deviate from the bulk value as the effective film thickness h ≲ hc ≈ 20σ related to the intrinsic properties of polymer films, while the Tg reduction is stronger for the free-standing films. However, a detailed study of Tg in the layers of polymer films, the distribution of entanglements inside the films, and the mobility of polymer chains near the surface is needed for the further understanding of the Tg reduction.
SUPPLEMENTARY MATERIAL
See the supplementary material for the simulation models, the explicit formulas of theoretical predictions based on Silberberg’s and the extended Silberberg models, and the determination of Tg from the total potential energy Utot(T) change.
We are grateful to Denis Andrienko for a critical reading of the manuscript. This work was supported by the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ERC Grant No. 340906-MOLPROCOMP. We also gratefully acknowledge the computing time granted by the John von Neumann Institute for Computing (NIC) and provided on the supercomputer JUWELS at the Jülich Supercomputing Centre (JSC) and the Max Planck Computing and Data Facility (MPCDF).
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
H.-P.H. performed the molecular dynamics simulations, analyzed the data, and wrote the paper. K.K. analyzed the data and wrote the paper.
Hsiao-Ping Hsu: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (equal); Methodology (equal); Software (lead); Visualization (lead); Writing – original draft (equal); Writing – review & editing (equal). Kurt Kremer: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Methodology (equal); Resources (lead); Writing – original draft (equal); Writing – review & editing (equal).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.