I. INTRODUCTION
Ultra-thin polymer films with thicknesses in the range of tens of nanometers are increasing in relevance for a wide array of applications in the field of nanotechnology and bioengineering, such as nanoelectronics, Reference Ito1 nanocomposites, Reference Ajayan, Schadler and Braun2 coatings, Reference Bertrand, Jonas, Laschewsky and Legras3 polymeric membranes, Reference Rowe, Freeman and Paul4 and biosensors. Reference Chen, Miao, He, Wu and Li5 For instance, in nanoelectronics, the thermal and mechanical responses of nanopatterned polymeric structures obtained from lithography are important in determining the stability and reliability of these systems. Reference Yoshimoto, Stoykovich, Cao, de Pablo, Nealey and Drugan6,Reference Sharma and Reiter7 In thin-supported and free-standing films, the effects of the free surfaces on the molecular structure and dynamics of polymer chains fundamentally change the behavior of these systems compared to the bulk, and these effects become significantly more pronounced as the film thickness is reached below 100 nm. The free-surface effects are complex and depend on many factors such as molecular weight of polymers, thickness of film, and polymer chemistry, but a common observation is that the presence of surfaces in polymer thin films may result in a notable reduction of thermomechanical properties by shifting the relaxation time and enhancing the molecular mobility. Reference Roth and Dutcher8–Reference Hanakata, Douglas and Starr12 Understanding the fundamental mechanisms pertaining to the thermomechanical properties, in particular, the glass transition temperature (T g) and elastic modulus (E) of sub-100 nm polymer thin films are important for numerous technological applications of soft nanomaterials, and are being studied through novel experimental and computational techniques suited for the length scales of interest.
Previous studies have shown that intermolecular noncovalent weak interactions between polymer chains are crucial in determining the thermomechanical responses of amorphous polymers. Experimental and theoretical studies have consistently shown that intermolecular forces between polymer chains can be tailored to tune the T g of bulk amorphous polymer. Reference Ito1,Reference Xia and Keten13 Concerning the mechanical behaviors, the investigations of energy partitioning have shown that the elastic modulus and yield response of amorphous polymer are mainly dominated by intermolecular noncovalent interactions (i.e., van der Waals interaction), Reference Hossain, Tschopp, Ward, Bouvard, Wang and Horstemeyer14,Reference Li, Mulder, Vorselaars, Lyulin and Michels15 and these mechanical properties are intimately related to the segmental dynamics of chains and the physical aging phenomena. Reference Lyulin, Balabaev, Mazo and Michels16–Reference Hudzinskyy, Michels and Lyulin19 Low molecular weight plasticizing and antiplasticizing additives are often used to soften or stiffen the polymer by changing its free volume and molecular packing. Reference Landel and Nielsen20,Reference Riggleman, Douglas and de Pablo21 Experimental investigations and simulations have suggested that the interface and free-surface effects on the thermomechanical behaviors of polymer films can be influenced by numerous factors associated with intermolecular forces, such as retained solvents, Reference Kim, Mundra, Roth and Torkelson22,Reference Torres, Stafford and Vogt23 monomeric interactions, Reference Xia and Keten13,Reference Ellison, Mundra and Torkelson24,Reference Xia, Hsu and Keten25 and chain-tacticity-related density variations. Reference Grohens, Brogly, Labbe, David and Schultz26,Reference Grohens, Hamon, Reiter, Soldera and Holl27
More recent efforts have been engaged in understanding the molecular origins of free-surface effects on the thermomechanical properties of polymer films. The thickness-dependent behavior of glass transition temperature and elastic modulus of polymer films have been characterized via experiments and simulations. Below T g, the elastic moduli of polystyrene (PS) and poly(methyl methacrylate) (PMMA) films and nanostructures decrease dramatically below ∼40 nm measured by a buckling-based metrology. Reference Stafford, Vogt, Harrison, Julthongpiput and Huang28 Torres et al. show that the T g and elastic moduli of sub-100 nm methacrylate polymer films both decrease with film thickness, but the decrease in the modulus of the films is not correlated with the observed decrease in T g. Reference Torres, Stafford and Vogt11 A study by Declambre et al. shows that the elastic modulus of free-standing antiplasticized PMMA nanostructure decreases below the bulk value when the width is below 50 nm measured by the beam-bending deformation. Reference Delcambre, Riggleman, de Pablo and Nealey29 On the other hand, recent experimental evidences show surface stiffening as observed in polymer thin films in the rubbery regime above T g, Reference O'Connell and McKenna30–Reference O'Connell and McKenna32 in which the underlying physics of mechanical behaviors could be fundamentally different from that below T g.
However, due to the difficulty in measuring the mechanical responses of the polymer free-standing thin film with free surfaces only (no substrate) as illustrated in Fig. 1(a), current experimental techniques of measuring the elastic modulus of the film inevitably introduce substrate and retained solvents effects. Those effects have been recently found to be crucial in determining the surface dynamics at nanoscale. Reference Evans, Narayanan, Jiang and Torkelson33 More importantly, how intermolecular forces between polymer chains interplay with the free surface that changes influence in T g and mechanical responses of polymer films remains challenging to discern, and the physics of this intriguing shift is yet to be fully explored. Nevertheless, a direct comparison of polymers with different chemistry via experiments and all-atomistic simulations is difficult since these parameters cannot be continuously and independently varied. Simplified potentials in coarse-grained molecular dynamics (CGMD) simulations of free-standing film with free surfaces provide a clear benefit in this regard in addition to computational efficiency.
(a) Illustration of polymer nanostructures with a width in nanoscale
$\left( {W \ll L,H} \right)$
(left) and the configuration of all-atomistic polymer free-standing film in simulations with periodic boundary conditions in x and y axes. (b) The mapping from the all-atom model of poly(methyl methacrylate) (PMMA) to the coarse-grained model (left) and the snapshot of the coarse-grained free-standing polymer film with two free surfaces on the top and bottom faces (right).
To facilitate a deeper understanding of the interplay between the intermolecular forces and free-surface effects on changing the thermomechanical responses of polymer thin films, we employ CGMD simulations to investigate independent effects of these parameters. For this purpose, we establish a CG approach that builds on realistic parameters derived from an all-atomistic MD simulation of PMMA. By employing this CG model as a starting point, we carry out large-scale MD simulations to assess the free-surface effects by examining the thickness-dependence behaviors of T g and modulus of the free-standing thin films with two free surfaces. For the modulus of the polymers, we focus only on the elastic region at a glassy state below T g, where polymers are deformed significantly faster than the intrinsic relaxation time, resulting in mechanical properties that are almost independent of temperature. The interplay between intermolecular forces and free-surface effect on film T g and modulus is evaluated by varying the depth of the effective weak interaction potential well between monomers. The interfacial and characteristic length scales that correlate with the scales where T g and E of the thin films deviate from the bulk value are examined and compared. A bilayer composite model derived from the spatial stress distribution is used to explain the physical origin of the size effect on the film modulus.
II. METHODS
A. Coarse-grained model development
To simulate films with a characteristic dimension (i.e., thickness) reaching ∼100 nm, we first present a coarse-grained (CG) modeling approach based on the output of all-atomistic simulations of PMMA. The one bead per repeat unit mapping from the all-atomistic model to the CG model and the resultant CG free-standing film configuration are illustrated in Fig. 1(b). The bonded interaction terms including the bond stretch and angle bending terms of the CG model are derived to match the corresponding probability distribution functions of the all-atom model. For the intermolecular noncovalent interactions between chains (including interchain and intrachain interactions), we employ the truncated Lennard-Jones (LJ) 12–6 potential:
where ε is the depth of the potential well, and σ corresponds to the radial distance where the potential energy is zero. These two parameters are calibrated to match the experimental density at 300 K and T g of the bulk PMMA. The resultant values of σ and ε are 6.5 Å and 0.3 kcal/mol, respectively, which yield a density of 1.17 g/cm3 at 300 K and a bulk T g of 388.6 K for our CG model. The cutoff distance of the intermolecular interaction R c is chosen to be 15 Å. A detailed description of the coarse-graining strategy and simulation protocols can be found in our previous work. Reference Xia, Mishra and Keten34
Utilizing the atomistically informed CG model as the generic polymer model, we are able to explore the effect of intermolecular noncovalent interactions between chains on the thermomechanical properties (i.e., T g and elastic modulus) of film systems with the existence of free surfaces at nanoscale. Taking the calibrated force–field parameters as a reference, we systematically vary ε from 0.2 to 0.4 kcal/mol to generate polymer models with different intermolecular interaction strength while keeping σ the same and thus the density nearly constant for all simulations. It should be noted that once ε is varied, it is no longer a representative of the original PMMA CG model. The similar two-bead type of CG models of methacrylate polymers has been recently reported. Reference Hsu, Xia, Arturo and Keten35 Comparing with the two-bead models, the bonded interactions in our one-bead model preserve the general backbone features of methacrylate polymers, and the different values of ε can be a representative of different interactions resulting from the different types of side-chain groups. The simpler single-bead model with reduced degrees of freedom and atomic friction gives rise to overall softer mechanical responses than more realistic two-bead CG models and all-atomistic models, Reference Rosch, Brennan, Izvekov and Andzelm36 hence the study here is meant to be largely comparative and qualitative.
B. Polymer free-standing thin film configuration
CGMD free-standing thin film simulations are carried out using 50 to 500 chains with 100 repeat units per chain, resulting in thin films with thicknesses ranging from ∼5 to ∼90 nm (along the z axis). Periodic boundary conditions of the simulation box are applied in the x–y plane with dimensions of 9 × 9 nm except for ∼5 nm thick films, in which we use larger dimensions of 18 × 18 nm in the x–y plane to improve the sampling. The film thickness is estimated as the difference between the maximum and minimum z-coordinates of the film. To create two free surfaces of the film, the box dimension along the z-axis is taken to be larger than the film thickness and the boundary condition in that direction is nonperiodic. The initial film structures are built based on a self-avoiding random walk algorithm and placed at the center of the simulation box along the z-dimension. The film first undergoes two annealing cycles from 550 to 210 K, and is then equilibrated at 250 K for 2 ns under the isothermal–isobaric ensemble (NPT). The time step is chosen to be 4 fs.
C. Glass transition calculation
The calculations of the glass transition temperatures of both the bulk and film samples are carried out using the method described by Tsige and Taylor Reference Tsige and Taylor37 from the measurement of the mean-squared displacement (MSD) data of monomers in polymer chains. The MSD of polymer beads is evaluated as:
where r
i(t) is the position of the ith bead at time t and
$\left\langle \ldots \right\rangle$
denotes the ensemble averages. The time-averaged values of g
0 from 40 to 2000 ps are plotted against temperature from 490 to 190 K at an interval of 20 K, which provides a metric of segmental mobility as a function of the temperature. The MSD data in the lower temperature region and higher temperature region can be fitted linearly, and the intersection of the two fits marks the T
g.
D. Calculation of elastic modulus
To calculate the elastic modulus, a strain-controlled uniaxial tensile deformation is simulated in the x direction (along the plane of the film), with a constant strain (engineering strain) rate of 0.5 ns−1, while the pressure in the transverse direction (y) is controlled to be zero using the Nose–Hoover barostat. The z dimension of the simulation box is automatically adjusted to encompass all the beads during the simulations. The stress component in the tensile direction is calculated based on the atomic virial stress tensor Reference Allen and Tildesley38 :
$$\tau _{{\rm{ij}}} = - {1 \over V}\left[ {\mathop\sum\limits_\alpha {m_\alpha \left( {v_\alpha } \right)_{\rm{i}} \left( {v_\alpha } \right)_{\rm{j}} \,- \mathop\sum\limits_{\alpha > \beta } {{{\partial U} \over {\partial r_{\alpha \beta } }}{{\left( {r_{\alpha \beta } } \right)_{\rm{i}} \left( {r_{\alpha \beta } } \right)_{\rm{j}} } \over {r_{\alpha \beta } }}} } } \right]\quad ,$$
where V is the volume of the actual polymer system instead of the simulation box, and m α and v α denote the mass and velocity of αth monomer, respectively. The distance between the monomer pair α and β is denoted by r αβ. U is the total potential energy including bond, angle, and pair-wise nonbonded interactions. While the strain rate used in our study is higher than that used in experimental studies, it is reasonably low for MD simulations, and the glassy modulus dependence on the strain rate is relatively weaker than rubbery modulus measurements. Reference Mulliken and Boyce39 Thus, qualitatively, we do not anticipate strong rate dependence effects on the comparative analyses presented herein. The elastic modulus is measured as the slope from 0 to 4% strain over 10 replicas.
III. RESULTS AND DISCUSSION
We first characterize how T g and elastic modulus E depend on the intermolecular interaction strength (ε) for bulk polymer systems (no free surfaces). The simulation results of the bulk T g and E are summarized in Table I. Both T g and E systematically increase with ε. Recently, we have shown that increasing ε from 0.2 (soft) to 0.4 kcal/mol (hard) leads to a monotonic increase in bulk T g, Reference Xia and Keten13 which can be ascribed to the fact that a stronger intermolecular interaction leads to a greater potential energy barrier for polymer chains to achieve cooperative segmental motion. The result of uniaxial tensile tests for glassy bulk polymer at T = 250 K is shown in Fig. 2. From the stress–strain responses, there is a distinct difference of deformation mechanisms between different ε values. For stronger intermolecular interaction (ε = 0.4 kcal/mol), four regimes can be observed: elastic, yield, strain softening, and hardening regions. This is similar to the experimental stress–strain observations on glassy polymers. Reference Mulliken and Boyce39 However, for weaker intermolecular interaction (ε = 0.2 kcal/mol), there is no obvious strain softening or hardening regime. The strain softening and hardening regimes with a stronger ε can be attributed to an increased activation energy associated with the molecular arrangement and interchain sliding after elastic deformation. Reference Rottler18 For the elastic response, the elastic modulus of the bulk polymer E increases with ε, which is consistent with previous findings, since the elastic work is primarily stored as nonbonded internal energy. Reference Hossain, Tschopp, Ward, Bouvard, Wang and Horstemeyer14
Comparison of glass-transition temperature (T g) and elastic modulus (E) of bulk polymer at a glassy state (T = 250 K) for different intermolecular interaction strengths (ε).

Stress–strain response for bulk polymers with different intermolecular interactions (ε = 0.2, 0.3, and 0.4 kcal/mol) under uniaxial tension at 250 K. The elastic modulus is determined by fitting the data with the slope from 0 to 4% strain over 10 replicates.
In polymer thin films, the regions near the free surfaces are characterized by an enhanced mobility and a reduced density compared to the interior bulk-like regions due to the less constraint on the surfaces. These regions are often called interfacial layers. The thickness of the interfacial layer, or the so-called interfacial thickness h
int, of a free-standing film is calculated on the basis of the Gibbs dividing surface (GDS) concept in thermodynamics, which can be used as an indication of the free-surface effect on the mobility of polymer films. Figure 3(a) shows the density profile along the z position of the film. The value of h
int is calculated by the equation:
$h^{{\rm{int}}} = {{\left( {h_{{\rm{max}}} - h_{{\rm{eff}}} } \right)} \mathord{\left/ {\vphantom {{\left( {h_{{\rm{max}}} - h_{{\rm{eff}}} } \right)} 2}} \right. \kern-\nulldelimiterspace} 2}$
, where h
max is the maximum film thickness defined by the distance between the maximum and minimum z coordinates of the film and h
eff is the effective thickness between the two GDS positions as labeled by the vertical dotted lines. Figure 3(b) shows the result of h
int as a function of film thickness for different ε at 250 K (below bulk T
g). From the result, it can be clearly observed that the polymer film with a lower value of ε has a larger h
int and thus a greater free-surface effect on the chain mobility for all the film thicknesses. The obtained interfacial thickness is about 0.5–1.0 nm, which is in a range similar to that reported by earlier MD studies of polymer films.
Reference Xia and Keten13,Reference Li and Strachan40,Reference Huang, Yang, An, Li and Chan41
For each ε, the film with a different thickness exhibits almost the same value of h
int. By taking the mean value of h
int over different film thicknesses, we obtain the average interfacial thickness for each ε, as listed in Table II.
(a) The density profile along the height of free-standing thin film (z) for ε = 0.3 kcal/mol at 250 K. (b) Interfacial thickness of mobile free-surface layer h int as a function of film thickness. The dashed lines indicate the mean value of h int over different film thicknesses.
Summary of different resultant surface-related properties that are used to quantify the free-surface effects and their relevant length scales as described in the text for different intermolecular interaction strengths (ε).

The thickness dependence of
$T_{\rm{g}}^{{\rm{film}}}$
of free-standing films with different intermolecular interactions is shown in Fig. 4. The
$T_{\rm{g}}^{{\rm{film}}}$
increases with the film thickness and then converges to their bulk Tg
for each ε. Over different thicknesses, the film with a larger value of ε exhibits a greater
$T_{\rm{g}}^{{\rm{film}}}$
. However, for a polymer with a larger ε, the
$T_{\rm{g}}^{{\rm{film}}}$
converges to its bulk value faster with increasing thickness compared to a softer polymer with a lower ε. The free-surface-induced reduction in the film T
g, which can be attributed to the enhanced mobile layer at the free surface as discussed in earlier studies,
Reference Sharp, Teichroeb and Forrest42–Reference Keddie, Jones and Cory45
is more significant for weaker intermolecular interactions. Therefore, the simulation result suggests that the free-surface effect on the depression of film T
g is stronger for weaker intermolecular interactions.
The thickness dependence of the glass transition temperatures of free-standing thin film
$T_{\rm{g}}^{{\rm{film}}}$
for different ε. The solid line is calculated as the best fit of Eq. (4), which results in the fit parameter of the characteristic thickness δT for
$T_{\rm{g}}^{{\rm{film}}}$
.

To further quantify the film thickness effects on T
g scaling and its relation to intermolecular interactions as shown in Fig. 4, we use the empirical fitting formula introduced by Kim et al.
Reference Kim, Jang and Zin46
to describe
$T_{\rm{g}}^{{\rm{film}}}$
:
where h is the film thickness,
$T_{\rm{g}}^{{\rm{bulk}}}$
is obtained from the bulk simulations (Table II), and δT is called the characteristic thickness that determines the saturation rate of
$T_{\rm{g}}^{{\rm{film}}}$
growth. The solid lines in the plot represent the best fit of Eq. (4) to the data, and the obtained result of δT is summarized in Table II. By employing this simple empirical form that involves only one fitting parameter, we are able to quantitatively assess how the relative free-surface effects on the film change with different intermolecular forces. From Table II, it can be observed that the value of δT depends on the intermolecular interaction ε and decreases with increasing ε. This trend is qualitatively consistent with h
int obtained from the thermodynamic analysis of the film, even though the value of h
int is smaller than δT for each ε. This result indicates that the thermodynamic analysis of free surface (GDS) can be well correlated with the thickness dependence of the film glass transition behavior and suggests that both h
int and δT relate to the depth of the free-surface layer with enhanced mobility for the systems studied herein.
Next, we measure the elastic modulus of free-standing films E film with different ε, as well as different thickness values. Figure 5(a) shows the snapshots of MD simulations during the tensile tests of the free-standing film. The thickness dependence of the modulus of the films E film at 250 K is shown in Fig. 5(b). Similar to the glass transition behaviors of the films, the free-surface-induced decrease in E film compared to the bulk value becomes more pronounced as the thickness decreases.
(a) MD simulation snapshots during the tensile deformation of the free-standing film with a thickness of 18 nm. The arrows indicate the direction of deformation of the film. (b) The elastic modulus of free-standing film E
film as a function of film thickness for different ε. The dotted line is the fit of empirical scaling formula in Eq. (5). The solid line is the prediction of Eq. (7) derived from the bilayer composite model. (c) Comparison of the empirical scaling results for the normalized T
g
$\left( {{{T_{\rm{g}}^{{\rm{film}}} } \mathord{\left/ {\vphantom {{T_{\rm{g}}^{{\rm{film}}} } {T_{\rm{g}}^{{\rm{bulk}}} }}} \right. \kern-\nulldelimiterspace} {T_{\rm{g}}^{{\rm{bulk}}} }}} \right)$
and E
$\left( {{{E^{{\rm{film}}} } \mathord{\left/ {\vphantom {{E^{{\rm{film}}} } {E^{{\rm{bulk}}} }}} \right. \kern-\nulldelimiterspace} {E^{{\rm{bulk}}} }}} \right)$
as a function of film thickness. (Inset) Comparison of normalized E for different ε.
It appears from our results that both E
film and
$T_{\rm{g}}^{{\rm{film}}}$
are influenced by the presence of the free-surface effects in free-standing thin films. This finding leads to the following question that we seek to answer here: does the characteristic length scale at which the deviation of E
film emerges from the bulk value correlates with that for the glass transition temperature?
As a first attempt to elucidate this issue, here we directly compare the free-surface effects on the modulus and T
g of free-standing thin films by employing the functional form of Eq. (4) for both
$T_{\rm{g}}^{{\rm{film}}}$
and E
film. Specifically, we describe the thickness dependence of film modulus as:
where δE is the characteristic thickness obtained from the modulus measurement. The best fit of Eq. (5) is denoted by the dotted lines in Fig. 5(b).
In essence, δE and δT can be considered as the comparable figures of merit that quantify the dependence of the elastic modulus and glass-transition temperature on film thickness, which is directly related to the length scale where the film properties deviate from the bulk values. The simulated values of δE as a function of intermolecular interaction strength are listed in Table II. The value of δE increases from 4 to 10 nm as ε decreases from 0.4 to 0.2 kcal/mol. This is reminiscent of the observation that the free-surface effect on T
g is greater for lower ε. The values of δE are greater than those values of δT for the same ε, which indicates that the convergence of the modulus to the bulk value with an increase in the film thickness occurs more slowly compared with
$T_{\rm{g}}^{{\rm{film}}}$
as shown in Fig. 5(c).
Our observations over a broad parameter range hint toward a general applicability of some of the concepts arising from previous experimental studies. A study by Torres et al. Reference Torres, Stafford and Vogt11 has shown that polymethacrylate films exhibit a much larger length scale for observed changes in the modulus at room temperature in comparison with the apparent T g. Forrest et al. Reference Forrest, Dalnoki-Veress and Dutcher47 also showed that the high-frequency mechanical properties of free-standing PS films measured by Brillouin light scattering (BLS) are not correlated with the film T g as the film thickness is varied. Wang and Mckenna Reference Wang and McKenna31 report that the significant T g reduction for PS when the film thickness is less than 20 nm using the liquid dewetting method. Additionally, earlier studies by McKenna's group and others Reference O'Connell and McKenna32,Reference Bodiguel and Fretigny48,Reference Hutcheson and McKenna49 have shown that the observed changes in rheological properties of polymer films with varying film thicknesses could not be well correlated with the T g measurement.
To further interpret our results and quantify the length scale of perturbation of the free surface, the information of distributions of the local stress and relaxation associated with the segmental mobility of polymer chains within the film can provide an insight into the nature of the surface effect. Figure 6(a) shows the plot of both spatial distribution of local stress and segmental relaxation time as a function of film position z. The local stress
$\tau (z) = {1 \over 2}\left[ {\tau _{xx} (z) + \tau _{yy} (z)} \right]$
is measured from the center to the free surface of the film intermediately after applying a biaxial step strain of 10% in x and y directions for ε = 0.2 kcal/mol. It can be clearly seen that the unrelaxed local stress is zero at the free-surface boundary, and increases more or less linearly with the depth into the film, as labeled by the solid slope in the plot. At ∼7 nm away from the free surface, the local stress saturates to the value in the interior region, which is indicated by the horizontal dashed line. Since the stress developed at a constant strain is related to the local modulus of the film, this analysis reveals that the free-surface effect that gives rise to a linear gradient in the modulus starting from the surface can penetrate several nanometers into the film.
(a) Spatial distributions of the stress and relaxation time across the film with a thickness of ∼78 nm from the center to the free surface (ε = 0.2 kcal/mol). The local stress
$\tau = {1 \over 2}\left( {\tau _{xx} + \tau _{yy} } \right)$
is measured in a biaxial deformation with a step strain of 10% in x and y directions. The local relaxation time is determined based on Eq. (6). The solid slope and dashed vertical lines correspond to the linear gradient of the stress within the surface layer and the average stress and relaxation time in the interior region of the film, respectively. h
surf is defined as the thickness of the linear slope region of the stress at the free surface labeled by the colored shadow in the plot. (b) The schematic of the bilayer composite model that consists of two soft surface layers with a thickness of h
surf and an interior bulk-like layer derived from the local stress distribution.

The local relaxation time τ0 plotted in Fig. 6(a) can be defined as the time needed by the MSD of CG beads within the local region to reach the bead size σ = 6.5 Å:
Compared to the stress distribution, it can be observed that the length scale of the relaxation gradient within the surface region is slightly shorter than that of the stress gradient. Considering the fact that the film elastic modulus converges to the bulk value at a larger film thickness compared to T g, this finding may suggest that length scales of film T g and modulus may not necessarily be the same because the length scales of the gradient of the local stress and mobility in the surface region can be different. It should be noted that it is very challenging to accurately predict the spatial T g within very small regions less than ∼1 nm using the current MSD approach. This is mainly because the segmental mobility of polymer chains can be very high and the monomers can diffuse from their original spatial region into other layers at a temperature close to or above T g. However, the local mobility via the measurement of local MSD should be related to T g, since both properties are measured based on the information of segmental dynamics in our study.
On the basis of the local stress distribution across the film, we employ a composite bilayer model as illustrated in Fig. 6(b) to explain the thickness dependence of E
film. The film with a total thickness of h can be considered as a composite with two softer surface layers with a thickness of h
surf and one stiffer bulk-like interior region based on the spatial stress distribution. The shape of the profile in spatial stress is similar to that in density, but the value of h
surf is larger than the interfacial thickness h
int determined from GDS. Based on the observation that h
int is nearly independent of film thickness, one can assume that this linear gradient holds for different film thicknesses. Following this, the surface layer has an average modulus of
$E^{{\rm{surf}}} = E^{{\rm{bulk}}} /2$
. As the distance from the free surface increases beyond the linear gradient region, the free-surface effect vanishes, and the modulus can be considered to be equal to E
bulk. If the film thickness h is less than 2h
surf, this means that the effects of the two free surfaces are pervasive throughout the film and there is no interior bulk-like region. This picture is qualitatively similar to the concept of the glassy bridge-reinforced model
Reference Merabia, Sotta and Long50–Reference Papon, Merabia, Guy, Lequeux, Montes, Sotta and Long52
that is introduced to explain the mechanical reinforcement of nanofilled elastomers. Then, E
surf becomes a linear function of h and is directly equal to E
film. Therefore, the E
film can be obtained based on the rule of mixtures:
$$E^{{\rm{film}}} = \left\{ {\matrix{ {E^{{\rm{bulk}}} \left( {1 - {{h^{{\rm{surf}}} } \over h}} \right)} & {h \ge 2h^{{\rm{surf}}} } \cr {{{E^{{\rm{bulk}}} } \over {4h^{{\rm{surf}}} }}h} & {h \lt 2h^{{\rm{surf}}} } \cr } } \right.\quad {\rm{.}}$$
The prediction of Eq. (7) is represented by the solid curves in Fig. 5(b). For thinner films with thicknesses below ∼10 nm, the prediction slightly underestimates the simulation results, which can be attributed to the lower accuracy in the assumption of linear modulus gradient in the surface region when the film thickness is less or comparable to 2h surf. However, without any input fitting parameters, the thickness dependence behavior of E film is well captured by the bilayer composite model. The values of h surf determined from the local stress distribution increase with decreasing ε (Table II), which is consistent with our previous empirical scaling analysis using Eq. (5). Even though the values of h surf arising from the stress perturbation of free surface are only in the range of 4–7 nm for different ε studied herein, the resultant change in modulus can be observed in a film with a thickness up to around 100 nm that is significantly larger compared to h surf. These values are also comparable with experimental observations. In the previous experimental study by Torres et al., Reference Torres, Stafford and Vogt11 the thickness of the free-surface layer for PMMA films on PDMS is reported to be about 3.5 nm. Our simulation prediction of h surf for ε = 0.4 kcal/mol is quite close to their value. In a recent study by Brinson et al., Reference Watcharotone, Wood, Friedrich, Chen, Qiao, Putz and Brinson53 they show that the interphase in the small local region near the substrate with a higher local modulus results in an altered elastic modulus of the film up to hundreds of nanometers, which supports our observation. In the simulations by Yoshimoto et al., Reference Yoshimoto, Jain, Nealey and de Pablo54 a free-standing polymer film exhibits local softening regions with lower storage modulus at two free surfaces via the measurement of local dynamic mechanical properties, and the spatial modulus exhibits a nearly linear gradient in the surface region, which is very similar to what we observe herein. Therefore, our finding physically explains why the E film approaches the bulk value as the film thickness increases by employing the bilayer composite model, and highlights the importance of determining the length scale of the interphase region arising from the free surfaces in polymer films.
IV. CONCLUSION
In this study, we provide a comprehensive CGMD study that quantifies how the glass-transition temperature and elastic modulus at the glassy state of free-standing polymer thin films depend on the intermolecular interactions and film thickness. Both bulk and film systems exhibit a greater value of T g and modulus with a similar backbone structure but stronger intermolecular interaction. However, the free surface has a stronger effect on the T g and modulus for the film with a weaker intermolecular interaction at nanoscale. By comparing the results of T g and modulus of the film, we find that the convergence of the modulus to the bulk value with an increase in the film thickness occurs more slowly compared with T g. This longer-range free-surface effect on the modulus other than T g may be further supported by the observation of a greater length scale of the local stress distribution in the surface region than that of the local relaxation time. The observed decrease in the modulus of the film exists up to the thickness of over ∼100 nm. However, the measurement of local stress distribution induced by a step strain deformation reveals that there exists a soft region with a length scale on the order of several nanometers, which is much smaller than the critical film thickness where the film modulus converges to the bulk value. The thickness dependence of the film modulus is well captured by employing the concept of the bilayer composite model, which physically explains why the modulus increases with the film thickness. Our findings establish a comprehensive and analytical description of free-standing film elastic properties, and highlight the importance of determining the length scale of the interphase region arising from the free surfaces in polymer thin films.
ACKNOWLEDGMENTS
The authors acknowledge funding by the Army Research Office (award #W911NF-13-1-0241) and Department of Civil & Environmental Engineering and Department of Mechanical Engineering at Northwestern University. A supercomputing grant from Quest HPC System at Northwestern University is acknowledged.





