Glass transition temperature of (ultra-)thin polymer films

In a majority of applications (commodity), amorphous polymeric materials are in the glassy state. Because of that, the glass transition temperature T_g region is of crucial importance,¹ 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 T_g, 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.⁹ Similarly, Xu et al.¹⁰ 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 T_g of polymers are by differential scanning calorimetry (DSC)¹² or by thermomechanical analysis (TMA),¹³ which can lead to deviating results.⁹ 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 T_g 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 T_g 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 chains³⁸ 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.⁴¹ 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 N_e = 28 beads at the standard melt density of ρ₀ = 0.85σ⁻³. The reference temperature for equilibration and cooling is T = 1ɛ/k_B. Throughout the whole paper, the Lennard-Jones (LJ) units of length (σ), energy (ɛ), and time (τ) are used.

For equilibration and structural properties at T = 1ɛ/k_B, 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,¹¹ with attractive non-bonded interactions, allowing for free surfaces where the pressure is kept at^52–57, P = 0.0ɛ/σ³ 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ɛ/k_B and ρ = 0.85σ⁻³, 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 T_g^11,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 potential^41,59,60 are introduced. This is sufficient to investigate the generic conformational properties of the films.

The ESPResSo++ package^61,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 n_c chains of N monomers (n_cN = 10⁵ for N ≤ 100 and n_c = 1000 for N = 500 and N = 2000), ranging from unentangled (N < N_e) to highly entangled (N ≫ N_e) films. The film thicknesses h range from thick $(h>Rg(0))$ to very thin $(h<Rg(0))$ films, $Rg(0)$ being the root mean square (rms) radius of gyration for bulk melt chains. The smallest h ≈ 9.0σ < 2d_T (see Fig. 1), d_T = 5.02σ being the tube diameter.⁴³ 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.⁴² 

We first analyze the chain conformations for confined melts at T = 1.0ɛ/k_B 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 R_e,‖(N, h) and R_g,‖(N, h) remain unchanged almost down to $h≈2Re(0)$⁠, 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 $h≈Re(0)$⁠, this turns into a linear decrease with h for R_e,⊥(N, h) (R_g,⊥(N, h)), while weakly increases for R_e,‖(N, h) (R_g,‖(N, h)).³³ Earlier, in 1982, Silberberg³¹ 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 R_e,⊥(N, h) described by the scaling function $fs,⊥(x=h/Re(0)(N))$⁠; see the supplementary material; however, it does not affect the parallel component. This idealized picture works well in the thick film regime down to $h≈2Re(0)$⁠. Sussman et al.³⁵ 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 $fes,⊥(x=h/Re(0)(N))$⁠, 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, $Re(0)/Rg(0)≈6$, while the distribution of $Rg(0)$ is not exactly Gaussian anymore. Therefore, our data on R_g,⊥(N, h) slightly deviate from $fes,⊥(x=h/Re(0)(N))$⁠. Furthermore, a detailed analysis of our data in the parallel direction reveals that the scaling behavior of R_e,‖(N, h) and R_g,‖(N, h) is better described by f_‖(x) = 1 + c_hx² instead of f_‖(x) = 1 + c_hx, with $x=(h/ξ(0)(N))1/2$ as in Ref. 64. Here, $ξ(0)(N)=Rg(0)(N)/N̄≈const$ (see Fig. S1 in the supplementary material) is the excluded volume screening length. $N̄=ρ0(Rg(0)(N))3/N≈0.2563N1/2$ is the degree of interdigitation of different polymers, also called generalized polymerization index⁶⁴ in bulk melts.

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 analysis⁶⁵ 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 N_kink 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 N_kink 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.6N_e). However, both datasets deviate from the bulk below and around h = h_c ≈ 20.0σ. Obviously, h_c is not related to $Rg(0)(N)$ $(Re(0)(N))$⁠, 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.³⁶ 

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ɛ/k_B. 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ɛ/σ³. 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 P_zz 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,⁵⁴ which results in a cooling rate Γ = ΔT/Δt = 8.3 × 10⁻⁷ɛ/(k_Bτ). The temperature is reduced in steps of ΔT = 0.025ɛ/k_B from T = 1.0ɛ/k_B to 0.2ɛ/k_B with a relaxation time between each step of Δt = 30 000τ ≈ 13τ_e (⁠$τe=τ0Ne2$ being the entanglement time defined by the tube model⁷⁰ and reptation theory,^71,72 and τ_e ≈ 2266τ estimated from simulations^42,43), i.e., subchains of the order of N_e can relax easily at higher temperatures close to T = 1.0ɛ/k_B.

To determine T_g, we perform a hyperbolic fit⁷⁴ to the density ρ(T) change with temperature $ρ(T)=c−a(T−T0)−b/2(T−T0+(T−T0)2+4ef)$, where c, T₀, a, b, and f are the fitting parameters. The density is defined by ρ(T) = n_cN/(h(T)L(T))², 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, T_g is defined by either T_g = T₀ 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 $Tg(0)(N)$ of the corresponding bulk melts^11,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 $Tg(0)(N)$ decreases with N for N < 2N_e, while $Tg(0)(N)≈Tg∞=0.6718(44)$ for N ≳ 2N_e.

Figures 4(c) and 4(d) show the results of T_g(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, T_g(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 by^30, $Tg(h)=Tg(0)(1−(h0/h)δ)$ with an exponent δ = 2.0(1), consistent with δ = 1.8(2) for the supported PS films.^30,75 The characteristic length h₀ 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, h₀ also depends on the chemical composition of polymer chains.⁷⁶ The deviation of T_g(N, h) from $Tg(0)(N)$ 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 T_g(N, h) seems to be almost independent of chains’ length as observed from experimental studies,⁷⁶ unlike the value of T_g(N) in bulk. Alternatively, one can also determine T_g from the bilinear fits of the total potential energy U_tot(T) change^54,77 (see the supplementary material). The estimates of T_g for N = 2000 are shown in Fig. 4. The T_g reduction remains the same; however, the uncertainty is much larger.

In summary, we confirm that the dependence of T_g 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 T_g, S, and N_kink (for films of highly entangled polymer melts) start to deviate from the bulk value as the effective film thickness h ≲ h_c ≈ 20σ related to the intrinsic properties of polymer films, while the T_g reduction is stronger for the free-standing films. However, a detailed study of T_g 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 T_g reduction.