Molecular Dynamics Investigations of Polystyrene-Based Binary Thin Film Systems_Interfacial Properties and Mechanical Behavior by Iyandri_TilukWahyono

VIEWS: 7 PAGES: 114

									 Molecular Dynamics Investigations of Polystyrene-Based
  Binary Thin Film Systems: Interfacial Properties and
                 Mechanical Behavior


                             A Thesis


Presented in Partial Fulfillment of the Requirements for the Degree
   Master of Science in the Graduate School of The Ohio State
                            University

                                By

                 Coleman Alleman, B.A.(Honors)

          Graduate Program in Mechanical Engineering



                    The Ohio State University

                               2011




                Master’s Examination Committee:
                  Professor Somnath Ghosh, Advisor
                            L. James Lee
 c Copyright by

Coleman Alleman

     2011
                                     Abstract



   The implementation and development of novel nanoscale biomedical devices and

procedures depends critically upon the ability to manufacture affordable nanoscale

polymer structures. This, in turn, depends upon understanding of the properties

of the polymer materials at these scales, where often the material properties are

divergent from those observed at larger scales. In this thesis, a molecular dynamics

(MD) based approach is developed to analyze these properties at an atomic resolution.

Under the hypothesis that the emergence of new properties in nanoscale polymers is

due to the increased importance of interface effects, the primary focus of this work is

the differentiation of material behavior at various interfaces.

   The first study presented here analyzes a polystyrene (PS) - carbon dioxide (CO2 )

interface. The main result of this study is a quantification of the impact of CO2 on

the glass transition behavior of PS. It is shown that the introduction of CO2 depresses

the glass transition temperature Tg significantly and that this depression increases

with increasing CO2 pressure. For the highest CO2 pressure studied (7 MPa), the

observed Tg for the polymer is 320 K, more than 50 K below the value for bulk PS.

In the study of this binary system, a number of techniques are also developed for

the measurement and analysis of properties such as free volume and mobility. These

techniques are applied to quantify the dependence of various important properties of

the PS material on spatial location, temperature, and CO2 pressure.

                                          ii
   The results of the first study suggest that the bonding of PS should be facilitated

by the introduction of CO2 , possibly enabling bonding at near room temperature.

The second study presented here examines this possibility, using a MD model of

PS thin films to study the impact of CO2 on the structure of a symmetric PS-PS

interface during bonding at 300 K. The properties of the interface are used to analyze

the results of simulations which show that the strongest interface is produced for a

CO2 pressure of ∼ 2 MPa. It is determined that this strength is largely determined by

the development of chain segments that bridge the interface. The number of bridges

that develop is shown to be dependent on atomic mobility near the interface, which

is a maximum for PS in the presence of CO2 at ∼ 2 MPa.

   The final study presented in this work examines the effects of free surface and

silica substrate interface effects on the glass transition behavior of PS thin films.

The change in Tg for freestanding and silica-supported PS thin films is estimated

using an MD model. It is shown that a freestanding PS thin film roughly 20 nm in

thickness exhibits a depression in Tg of 33 K below the bulk value. It is also shown

that a silica-supported PS thin film exhibits a depressed Tg , but that the competition

between polymer-substrate and polymer free surface effects makes this depression a

pronounced 19 K for a film of identical composition. It is concluded that the deviation

in the thin film Tg values is driven by changes in atomic mobility due to interface

and free surface effects.




                                          iii
 Dedicated to my wife,

  Joanna Alleman,

for her love and support.




           iv
                              Acknowledgments



   I would like to begin by acknowledging my advisor, Professor Somnath Ghosh,

for his guidance and continued financial support throughout this work. His vision

for this work and his insights, conveyed over countless hours of discussion, have been

truly invaluable. My sincerest gratitude goes to Professor Ghosh. Without him, this

work would not have been possible.

   I would also like to acknowledge Professor L. James Lee, with whom I had many

fruitful conversations over the course of my research and whose experimental work

provided both the starting point for my research and the crucial validation for the

computational models. My thanks also to Professor Lee for serving as a member of

my examination committee.

   I would like to thank Prof. Sherwin Singer for the valuable insights I gained during

his course and our discussions thereafter.

   I take this opportunity to thank the past and current lab members and friends

at the Computational Mechanics Research Laboratory for the essential technical and

moral support they provided throughout the course of my study. Many thanks es-

pecially to Anand Srivastava, with whom I had inumerable discussions and shared

many hours of combined effort.

   This work has been supported by the NSF Division of Engineering Education and

Centers through grant no EEC-0425626 to the Center for Affordable Nanoengineering

                                             v
of Polymeric Biomedical Devices (CANPBD) at the Ohio State University (Professor

L. James Lee, Director). This sponsorship is gratefully acknowledged. Computing

support from the Ohio Supercomputer Center through grant PAS0191 is also grate-

fully acknowledged.




                                       vi
                                                             Vita



1981 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Born: Dayton, Ohio

2005 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.A. Philosophy,
                                                                                     The Ohio State University
2006-2008 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Student, Dept. of Mech. Engineering,
                                                                              The Ohio State University
2008-2010 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Graduate Research Assistant,
                                                                              The Ohio State University
2010-2011 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Graduate Fellow,
                                                                              The Ohio State University


                                                    Publications

Research Publications

Anand Srivastava, Coleman Alleman, Somnath Ghosh and L.J. Lee ”MD simula-
tion based evaluation of Tg of PS in the presence of carbon dioxide”, Modelling and
Simulation in Materials Science and Engineering, 18: 1363-1396, 2010.

Coleman Alleman, Anand Srivastava and Somnath Ghosh ”Understanding the Ef-
fects of Carbon Dioxide on Interfacial Bonding of Polystyrene Thin Films Through
Molecular Dynamics Simulations”, Journal of Polymer Science Part B: Polymer
Physics, In review.

Coleman Alleman, L. J. Lee and Somnath Ghosh, ”Silica-Polystyrene Interface
Effects on Local Tg ”, To be submitted.

Conference Publications



                                                                vii
”Nanoscale Modeling of CO2 -Assisted Bonding Processes in Polymer Thin Films,”
Coleman Alleman, Anand Srivastava and Somnath Ghosh, US Conference of The-
oretical and Applied Mechanics, June, 2010, University Park, PA”

”Computational Design Tools for Fabrication of Polymer Nanodevices,” Anand Sri-
vastava, Coleman Alleman and Somnath Ghosh, Ohio Innovation Summit, April,
2010, Columbus, OH”

”MD Simulation based estimation of glass transition temperatures of polystyrene in
the presence of CO2 ,” Anand Srivastava, Coleman Alleman and Somnath Ghosh,
US National Congress on Computational Mechanics, July, 2009, Columbus, OH.

”Multiscale Modeling Approaches for Nanotechnology,” Somnath Ghosh, Anand Sri-
vastava and Coleman Alleman, Invited talk at Nanocomputation symposium - Ohio
Innovation Summit, April, 2008, Dayton, OH.



                              Fields of Study

Major Field: Mechanical Engineering

   Studies in Molecular modeling : Prof. Somnath Ghosh




                                       viii
                                Table of Contents




                                                                                                             Page

Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .                          ii

Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .                         iv

Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .                             v

Vita . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .                         vii

List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .                         xii

List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .                         xiii

List of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii

1.   Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .                           1

      1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .                         .      1
      1.2 Methodology: Molecular Dynamics Modeling . . . . . . . . . . .                                 .      3
          1.2.1 Newtonian Mechanics and Interatomic Potential Functions                                  .      3
          1.2.2 Ensemble Averaging and Hamiltonian Mechanics . . . . .                                   .      4
          1.2.3 Simulation, Pre- and Post-Processing Algorithms . . . . .                                .      6
          1.2.4 System Representations . . . . . . . . . . . . . . . . . . .                             .      7

2.   Molecular Modeling: Polystyrene, Carbon Dioxide, and Silica                     . . . . . .                9

      2.1 Introduction . . . . . . . . . . . . . . . . . . .     .   .   .   .   .   .   .   .   .   .   .      9
      2.2 Potential Functions . . . . . . . . . . . . . . . .    .   .   .   .   .   .   .   .   .   .   .     10
          2.2.1 Potential Functions for Polystyrene . . .        .   .   .   .   .   .   .   .   .   .   .     10
          2.2.2 Potential Functions for Carbon Dioxide           .   .   .   .   .   .   .   .   .   .   .     14
          2.2.3 Potential Functions for Silica . . . . . .       .   .   .   .   .   .   .   .   .   .   .     17
      2.3 Modeling Component Interactions . . . . . . .          .   .   .   .   .   .   .   .   .   .   .     21

                                             ix
         2.3.1 Polystyrene - Carbon Dioxide Interactions . . . . . . . . . .            22
         2.3.2 Polystyrene - Silica Interactions . . . . . . . . . . . . . . . .        25
     2.4 Creation of Polystyrene Initial Configurations . . . . . . . . . . . .          26

3.   Initial Investigations: Interface Effects and Glass Transition in Finite
     Polystyrene-Carbon Dioxide Systems . . . . . . . . . . . . . . . . . . . .         28

     3.1 Inroduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .      28
     3.2 Experimental Methods For Determining Glass Transition Temperature              30
     3.3 System Parameters and Simulation Protocol . . . . . . . . . . . .              31
         3.3.1 Carbon Dioxide Bath . . . . . . . . . . . . . . . . . . . . .            33
     3.4 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . .       38
         3.4.1 Regions of Material Behavior: Bulk and Surface Properties .              39
         3.4.2 Effect of CO2 on Polystyrene Tg . . . . . . . . . . . . . . .             40
     3.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .      42

4.   Understanding the Effects of Carbon Dioxide on Interfacial Bonding of
     Polystyrene Thin Films . . . . . . . . . . . . . . . . . . . . . . . . . . . .     44

     4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .       44
     4.2 Adhesive Strength Development at a Symmetric Interface . . . . .               46
     4.3 Model Development and Simulation Procedures . . . . . . . . . . .              48
         4.3.1 Overview of Modeling Scheme . . . . . . . . . . . . . . . . .            48
         4.3.2 Model Materials and Configuration . . . . . . . . . . . . . .             49
     4.4 Investigation of Thin Film Behavior: Model Validation and Prelim-
         inary Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .      52
         4.4.1 Film Density Characteristics . . . . . . . . . . . . . . . . .           52
         4.4.2 Chain Dynamics in the Thin Films . . . . . . . . . . . . . .             53
     4.5 Numerical Experiments on the Bonding Process . . . . . . . . . . .             55
         4.5.1 Simulation of Bonding in PS-CO2 Binary Systems . . . . . .               55
         4.5.2 Density Profiles and Diffusion at the Interface . . . . . . . .            56
         4.5.3 Free Volume Characteristics and Chain Bridging at the In-
                 terface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .    59
         4.5.4 Dynamic Behavior of Chains at the Interface . . . . . . . .              62
         4.5.5 Simulation of Debonding . . . . . . . . . . . . . . . . . . . .          63
         4.5.6 Interface Integrity: A Strain Energy Density Based Metric .              64
     4.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .      66

5.   Quantification and Analysis of Effects of the Silica-Polystyrene Interface
     on The Glass Transition in A Supported Polystyrene Thin Film . . . . .             70

     5.1 Introduction     . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .   70


                                           x
      5.2 Simulation Details . . . . . . . . . . . . . . . . . . . . . . . . . . .     71
          5.2.1 Generation of Initial Configuration: Silica Slab . . . . . . .          71
          5.2.2 Generation of Initial Configuration: Polystyrene Thin Film .            73
          5.2.3 Creation of Binary Silica Substrate - Supported Polystyrene
                 Thin Film System . . . . . . . . . . . . . . . . . . . . . . .        74
          5.2.4 Simulation Protocol for Determination of Polystyrene Tg . .            74
      5.3 Simulation Results: Mobility-Based Polystyrene Tg Measurements .             75
      5.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .    79


Appendices                                                                             83

A.   Symplectic Integrators . . . . . . . . . . . . . . . . . . . . . . . . . . . .    83


B.   Determination of a Maximum Timestep for Stability of the Integration
     Algorithm in Molecular Dynamics Simulations . . . . . . . . . . . . . . .         85


Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .   87




                                            xi
                                List of Tables



Table                                                                                Page

2.1   Parameter values (energies in kcal/mol) in the potential energy func-
      tions for the all atom model of polystyrene. The C-H LJ parameter is
      calculated via the Lorentz-Berthelot mixing rule.[1] . . . . . . . . . .         13

2.2   Parameter values (energies in kcal/mol) in the potential energy func-
      tions for PS. The mixed LJ parameters are calculated via the Lorentz-
      Berthelot mixing rule.[1] . . . . . . . . . . . . . . . . . . . . . . . . .      15

2.3   Parameter values (energies in kcal/mol) in the potential energy func-
      tions for CO2 . The mixed LJ parameters are calculated via the Lorentz-
      Berthelot mixing rule. . . . . . . . . . . . . . . . . . . . . . . . . . .       16

2.4   Parameter values (energies in kcal/mol) in the potential energy func-
      tions for SiO2 . Note that there is only a coulombic part to the Si-Si
      interaction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .     19

2.5   Parameter values (energies in kcal/mol) in the TV potential energy
      functions for the silica slab. Here, OS refers to oxygen atoms in a
      silanol group. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .     20

2.6   Mixed LJ parameters for PS-CO2 interaction calculated via the CM-
      HHG mixing rule. Energies are in kcal/mol and distances are in
      angstroms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .     24

2.7   LJ parameters for PS-Silica interactions. Energies are in kcal/mol and
      distances are in angstroms. The atom types Os and Ox refer to atoms
      on the surface and within the silica slab, respectively . . . . . . . . .        25

3.1   Number of CO2 required to achieve given pressures at 300K . . . . . .            36



                                        xii
                               List of Figures



Figure                                                                                Page

2.1   Polystyrene monomer . . . . . . . . . . . . . . . . . . . . . . . . . . .         11

2.2   CO2 molecule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .        16

2.3   Representation of silica. . . . . . . . . . . . . . . . . . . . . . . . . .       17

2.4   A plot of the Si-O pair potential energy for the modified BKS potential. 18

2.5   Comparison of fv data from simulations of a binary polystyrene-CO2
      system with PALS results [2] at ambient pressure . . . . . . . . . . .            22

3.1   A schematic of the molecular dynamics setup for a finite polystyrene
      system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .     32

3.2   Calibration of CO2 pressure at 300K. . . . . . . . . . . . . . . . . . .          34

3.3   A schematic of molecular dynamics setup for finite PS-CO2 system
      with repuslive walls at the boundary. . . . . . . . . . . . . . . . . . .         35

3.4   Free volume fraction at 400 psi as a function of temperature for dif-
      ferent system sizes. Error bars represent the standard deviation of the
      measurement. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .        35

3.5   Evolution of density as the PS system is cooled from 348K to 323K for
      a duration of 1 ns (cooling rate of 0.05K/ps). The density converges
      after around 1 ns and further simulation does not cause any perceivable
      change in density. . . . . . . . . . . . . . . . . . . . . . . . . . . . .        37

3.6   Density variation across thickness at 400K and different pressures of
      CO2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .       38

                                         xiii
3.7   Gradients in density, free volume fraction and mean squared displace-
      ment for PS-CO2 systems at 400 psi. . . . . . . . . . . . . . . . . . .       39

3.8   In (a), specific volume is plotted against temperature for finite systems
      with different CO2 pressures. The glass transition point is identified by
      the change in slope (rate of thermal expansion). Standard deviation
      in measurements is shown by error bars. In (b), the determination of
      Tg from MD is compared with theory[3] and experiment.[4] The error
      in the MD measurements is bounded by the error shown in (a). . . .            41

4.1   Schematic representation of an experimental setup for the CO2 assisted
      bonding process in PS thin films for fabricating nanoscale devices. . .        45

4.2   A molecular simulation system for the thin film bonding process with
      CPK representation of CO2 : (a) The PS thin films are represented by
      dynamic bonds and are shaded differently for identification, (b) The
      surfaces of the PS thin films are represented by an interpolation of the
      van der Waals surfaces. . . . . . . . . . . . . . . . . . . . . . . . . .     48

4.3   Polystyrene density along the nonperiodic thickness direction for the
      15x200, 800psi thin film system. . . . . . . . . . . . . . . . . . . . . .     51

4.4   Polystyrene density along the nonperiodic thickness direction. A smooth
      curve is fitted to highlight the change in profiles upon addition of CO2 . 53

4.5   Percentage volume change for two PS thin film under various CO2 pres-
      sure conditions (zero pressure corresponds to absence of CO2 ). Error
      is of the order of the marker size. . . . . . . . . . . . . . . . . . . . .   54

4.6   Layer-wise mobility, MSDrel , for polymer atoms in a thin film. Error
      bars represent the standard deviation in the measurement. . . . . . .         55

4.7   Surface separation at the interface versus time for the 15 chain molecule,
      200 psi PS-CO2 system. The curve is typical of all simulation systems,
      except that the equilibrium separation S0 depends on CO2 pressure, as
      seen in Figure 4.8. . . . . . . . . . . . . . . . . . . . . . . . . . . . .   57

4.8   Equilibrium surface separation S0 at the interface following bonding
      simulations as a function of CO2 pressure. Error bars represent the
      standard deviation in the measurement. . . . . . . . . . . . . . . . . .      58


                                        xiv
4.9   Average monomer inter-penetration as a function of CO2 pressure for
      two different timescales and film thicknesses. Error bars represent the
      standard deviation in the measurement. . . . . . . . . . . . . . . . . .       59

4.10 Average number of cavities at the interface as a function of CO2 pres-
     sure for two different timescales and film thicknesses. Error bars rep-
     resent the standard deviation in the measurement (for the 33 ns sim-
     ulations, the error is smaller than the marker size). Lines are a guide
     to the eye. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .   60

4.11 Average volume of cavities at the interface as a function of CO2 pres-
     sure for two different timescales and film thicknesses. Error bars rep-
     resent the standard deviation in the measurement. Lines are a guide
     to the eye. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .   61

4.12 Number of chain segments bridging the interface as a function of CO2
     pressure for two different timescales and film thicknesses. Error bars
     represent the standard deviation in the measurement. Lines are a guide
     to the eye. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .   62

4.13 Mean-squared displacement over 200ps as a function of CO2 pressure.
     Error bars represent the standard deviation in the measurement. . .             63

4.14 Stress-strain curves for all debonding simulations. Carbon dioxide
     pressures are given in the legend. . . . . . . . . . . . . . . . . . . . .      64

4.15 Strain energy density developed up to failure as a function of CO2
     pressure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .   66

4.16 Strain energy density developed up to failure as a function of simulation
     time for the 10x200 polymer system under two CO2 pressure conditions. 67

5.1   Pictorial catalogue of SiO2 surface defects. In each figure, the darker
      atom(s) constitute the defect. The fictitious ’bonds’ are drawn to high-
      light the coordination of each atom. . . . . . . . . . . . . . . . . . . .     72

5.2   Representation of the creation of silanol groups on the silica slab sur-
      face. The dashed circle represents the original oxygen placement. . . .        73




                                        xv
5.3   Comparison of Tg calculations for bulk, freestanding thin film and
      supported thin film configurations. . . . . . . . . . . . . . . . . . . .        75

5.4   Measured values of Tg for the freestanding and supported films across
      a number of cooling rates. The dashed lines are fit to the data using
      Equation 5.2. The values of the coefficients are A=1.05±0.06 ps/K
      and B=181.0±12.6 K. The values of T0 are 342 K for the freestanding
      film, 356 K for the supported film, and 375 for the bulk configuration.           77

5.5   Comparison of mobility-based Tg calculations for the region near the
      polymer free surface and for the core region in a freestanding polystyrene
      film, and similar analyses for the free surface, interface, and core re-
      gions for a supported film of identical composition. . . . . . . . . . .        78

5.6   Layerwise mobility for freestanding and supported polystyrene films at
      two temperatures. The solid horizontal line corresponds to the MSD
      of the bulk system at the given temperature. For the supported film,
      the position of ∼ -11 nm corresponds to the surface of the silica slab
      substrate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .   79

5.7   Layerwise polymer densities for the freestanding and supported thin
      film systems. The horizontal solid line represents the density measured
      in the bulk system. . . . . . . . . . . . . . . . . . . . . . . . . . . . .    80

5.8   Layerwise Tg estimates for a 21 nm supported thin film. . . . . . . .           82




                                        xvi
                           LIST OF SYMBOLS



kB           Boltzmann constant = 1.381 x 10−23 [m2 kg.s−2 K−1 ]
MW           Molecular weight [g/mol]
Rg           Radius of gyration [˚]
                                 A
T            Thermodynamic temperature [K]
Tg           Glass transition temperature [Kelvin]
Tbg          Bulk Glass transition temperature [Kelvin]
Tig          Interface Glass transition temperature [Kelvin]
Tsg          Surface Glass transition temperature [Kelvin]
       1
β=    kB T
             Boltzmann factor [kj/mol]
ρ            Density [kg/m3 ]
σ            van der Waal radius of an atom [˚]
                                              A




                                      xvii
                     Chapter 1: INTRODUCTION




1.1    Motivation

   In the development of new medical techniques on the frontiers of disease detection

and treatment, nanoscale devices are increasingly playing a critical role. There are

numerous examples where such devices are sure to have a profound impact on patient

outcomes, including lab-on-a-chip devices for cancer detection and nanoparticle-based

systems for targeted drug delivery. Nanoscale polymer processing is critically impor-

tant for the fabrication of affordable biomedical nanodevices and remains a serious

challenge. The most serious limitation to fabrication is due to the increasing signifi-

cance of interface effects as the size of the devices approaches the nanoscale. At these

small size scales, the bulk properties of the constituent materials are less relevant,

and new material behavior emerges. Thus, it is impossible to design a device com-

ponent without considering the effects of the processing environment or interactions

with attached components.

   One of the most critical emergent phenomena in nanoscale polymer devices is

a new regime of temperature dependence for material properties such as stiffness.

The temperature threshold for the changeover in properties from disordered solid to

liquid-like behavior in polymers is known as the glass transition temperature (Tg ). It

                                          1
is well-established that there exists a surface layer at any polymer free surface, where

a depressed surface glass transition temperature (Ts ) can be distinguished from the
                                                   g

higher bulk glass transition temperature (Tb ). Additionally, some results suggest
                                           g

the existence of a differentiated layer at any interface, with an associated interfacial

glass transition temperature (Ti ), where the departure from Tb is governed by 1) the
                               g                              g

geometrical structure of the interface and 2) the interaction bewteen the polymer and

the material on the other side of the interface. This sort of differentiation between

different regions of material behavior is critically important for the manufacture of

polymer nanodevices for two reasons. First, as a polymer device becomes smaller, if

the region of interfacial properties remains roughly constant, the interfacial properties

are increasingly dominant in determining the properties of the entire device. Second,

nanoscale patterning or bonding can be more readily accomplished if the difference in

material properties in different regions can be exploited. For example, a softened free

surface may be beneficial for bonding two components together while a more solid

core prevents unwanted deformation in the rest of the part.

   From the above discussion, it is clear that there is a need for analytical techniques

that are able to probe material behavior at a very fine scale and that are able to differ-

entiate such behavior within a single material. Current experimental techniques are

limited to taking measurements of the behavior or to observing in a very limited way

the evolution of a non-symmetric material interface. A noted exception to this limi-

tation is a technique known as Positron Annihilation Lifetime Spectroscopy (PALS),

which can be used to probe polymer free volume in a fairly localizable way. However,

even PALS is not able to provide clear measurements of the localized deviations in

material behavior mentioned above. As such, molecular dynamics (MD) modeling is


                                           2
employed in the current study as a method of providing information about material

behavior at the atomic scale, which is sufficient to study localized phenomena like the

free surface induced depression in Tg .


1.2      Methodology: Molecular Dynamics Modeling

   In a MD model system, the positions and momenta (velocities) of individual atoms

are calculated at each timestep during a simulation. This section presents the com-

ponents of the models from which these quantities are calculated. A brief discussion

of the methods used to analyze the results of the MD simulations is also given.

1.2.1       Newtonian Mechanics and Interatomic Potential Func-
            tions

   The basis for the time evolution of a MD model system is Newton’s equation of

motion
                                                ∂p
                                          F =                                     (1.1)
                                                ∂t

relating force to change in momentum. This differential equation is solved in three di-

mensions for each timestep and for each atom in an n-atom ensemble. In the simplest

formulation of an MD simulation, this results in 6n first-order ordinary differential

equations to be solved for each timestep, one for each component of velocity and

position.

   In an equilibrium MD simulation, the force on an atom is governed by a (multi-

body) potential energy function as

                                                ∂U
                                      F =−                                        (1.2)
                                                ∂x




                                            3
        ∂U
where   ∂x
                 is the gradient of the potential energy function U. In turn, the potential

function U is defined for an atom i as


                                            U = Ub + Ua + Ut + Up                                           (1.3)


where Ub is the potential energy due to bond stretching, Ua is due to bond bending,

Ut is due to bond torsion, and Up is due to non-bonded pair interactions. In the

current study, the pair potential Up is assumed to be radially symmetric, so that the

force on an atom becomes
         N
                          ∂Ub        ∂Up                             ∂Ua                            ∂Ut
   F =              Bij        + Pij          +               Aijk         +               Tijkl            (1.4)
             j
                          ∂rij       ∂rij         j   k
                                                                     ∂θijk     j   k   l
                                                                                                   ∂φijkl

where the coefficients Bij , Pij , Aijk , and Tijkl are −1 if the interaction is defined

and 0 otherwise. The quantities rij , θijk , and φijkl refer to the bond length, bond

valence angle, and bond torsional angle, respectively. The functional forms of each

component of the total potential energy function are discussed in some detail in the

appropriate sections below.

1.2.2        Ensemble Averaging and Hamiltonian Mechanics

   The form of the governing equations for an MD system is complicated slightly

by the need to employ a particular ensemble sampling to obtain averages of various

thermodynamic quantities of interest. In a particular ensemble, two state variables

such as volume and temperature are held constant as the system evolves. In this

work two of the most common ensembles are employed. The first is the canonical

ensemble, where volume and temperature are held constant, and the second is the

isobaric-isothermal ensemble, where temperature and pressure are held constant. The

                                                                      e
development of the equations that follow is from a paper by Suichi Nos´. [5]

                                                          4
   To introduce the mechanisms by which the constant values are maintained, it is

convenient to recast the equations of motion in Hamiltonian form

                                  dqi   ∂H     pi
                                      =     =                                  (1.5a)
                                  dt    ∂pi   mi
                                  dpi    ∂H       ∂U
                                      =−      =−                              (1.5b)
                                  dt     ∂qi      ∂qi

where pi and qi are the momentum and position, respectively, of the i-th particle.

Here, the Hamiltonian H is given by

                                              p2
                                               i
                                 H=              + U(q)                         (1.6)
                                         i
                                             2mi

where mi is the mass of the i-th particle, and U(q) is the total potential energy for

the system. For the canonical ensemble, a time-scaling variable s is introduced, and

the Hamiltonian is extended to

                               p2
                                i
                   H=                + U(q) + p2 /2Q + 3NkB T ln(s)
                                               s                                (1.7)
                          i
                              2mi s2

Here, ps is the conjugate momentum for s, and Q is a corresponding ’mass’ for s. The

equations of motion for the extended system are

                       dqi    pi
                           =                                                   (1.8a)
                       dt    mi s2
                       dpi    ∂U
                           =−                                                 (1.8b)
                       dt     ∂qi
                        ds   ∂H     ps
                           =      =                                            (1.8c)
                        dt   ∂ps    Q
                       dps    ∂H   1               p2
                                                    i
                           =−    =                      − 3NkB T              (1.8d)
                        dt    ∂s   s          i
                                                  mi s2

For the isothermal-isobaric ensemble, a volume-scaling variable V is introduced, and

the Hamiltonian is extended to

                   p2i                                          p2
    H=                2/3 s2
                             + U(V q) + ps /2Q + 3NkB T ln(s) + V + Pex V
                                  1/3    2
                                                                                (1.9)
           i
               2mi V                                           2W

                                             5
where pV is the conjugate momentum for V , W is the corresponding ’mass’ term,

and Pex is the target pressure. This extended Hamiltonian gives rise to the following

equations of motion

                  dqi       pi
                      =                                                        (1.10a)
                  dt    mi V 2/3 s2
                  dpi    ∂U
                      =−                                                       (1.10b)
                  dt     ∂qi
                   ds   ∂H      ps
                      =      =                                                 (1.10c)
                   dt   ∂ps      Q
                  dps    ∂H   1                p2
                                                i
                      =−    =                          − 3NkB T                (1.10d)
                  dt     ∂s   s        i
                                           mi V 2/3 s2
                  dV   ∂H     pV
                     =      =                                                  (1.10e)
                  dt   ∂pV    W
                 dpV    ∂H       1              p2
                                                 i        ∂U
                     =−      =                   2/3 s2
                                                        −     qi   − Pex        (1.10f)
                  dt     ∂V    3V      i
                                            mi V          ∂qi

Integration of the equations of motion 1.8 or 1.10 with an algorithm that conserves the

associated Hamiltonian produces sets of coordinates that sample from the appropriate

ensemble, thereby allowing averages of various thermodynamic properties to be taken.

Such algorithms are known as symplectic integrators, and are briefly discussed in

Appendix A.

1.2.3    Simulation, Pre- and Post-Processing Algorithms

   To integrate the equations of motion for the various MD systems simulated in

this work, the programs DLPOLY and LAMMPS are used. [6, 7] Both porgrams em-

ploy efficient parallelization schemes; DLPOLY uses a domain decomposition method

to parallelize the computational load, and LAMMPS uses a spatial decomposition

method. In both programs, the second-order symplectic verlet algorithm is imple-

mented to integrate the equations of motion.


                                            6
   Pair potential functions are implemented within the MD program frameworks

as needed. Examples include the slow push-off potential for fast equilibration [8]

and the modified BKS and truncated Vessal potentials for the simulation of silica.[9]

Additionally, a number of pre- and post-processing algorithms are developed to create

atomic configurations, to analyze trajectories for atomic mobility, and to calculate Tg ,

to name just a few. The particular algorithms and their implementations are discussed

in detail in the appropriate sections.

1.2.4     System Representations

   Owing to the large number of equations and the fast dynamics of the typical

system, practical limits for MD simulations are for times on the order of nanoseconds

and for sizes on the order of nanometers. This corresponds to orders of 106 timesteps

for n ≡ 107 atoms. It can be seen, then, that MD is useful for studying the local

distributions of the types of material behavior discussed in the previous section only if

this behavior is localized to a small region and takes place on a short timescale. With

current computer technology, MD cannot be used to simulate a full device, or even

a device component of the type in question in full atomic detail. Thus, study of the

material behavior presented here is dependent on careful choices about representing

interesting phenomena in a limited material subset, such as that at the free surface

of a polymer.

   It is well known that the relaxation time for a realistic polymer is orders of mag-

nitude longer than that which can be simulated by MD methods. Phenomena such




                                           7
as the glass transition or interfacial bonding play out on timescales orders of magni-

tude longer again. In the present analyses, a great deal of attention is paid to these

limitations, and novel techniques are developed to circumvent these barriers.




                                          8
 Chapter 2: MOLECULAR MODELING: POLYSTYRENE,
                 CARBON DIOXIDE, AND SILICA




2.1    Introduction

   Polystyrene, carbon dioxide, and silica are ubiquitous materials, found in everyday

devices such as compact disc cases, soda water, and window glass. These materials

are also very important in the design and fabrication of nanoscale biomedical de-

vices. From the perspective of atomic structure, carbon dioxide and silica are among

the simplest materials to be found, and polystyrene, while slightly more complex, is

nonetheless composed of simple repeat units of 16 atoms each. The molecular mod-

eling of these simple forms of matter, on the other hand, is far from simple. The

formulation and execution of the potential functions governing the behavior of the

atomic constituents of these materials has been the subject of a considerable body of

work in quantuum chemistry and ab initio molecular dynamics.

   Here, the fruits of these labors are summarized - for each material and simulation

protocol, a particular potential field is selected from among those available. Modifica-

tions to the original potential functions are made to suit the needs of the simulations,

and these, too, are presented in this chapter. Additionally, investigation of the proper



                                           9
interactions between constituents is documented here, along with the final functional

forms employed in the various simulations.


2.2     Potential Functions

   In the molecular dynamics simulations of this study, the equations of motion

for the simulation systems are determined by calculating the forces exerted on each

atom by the other atoms in the system. These forces, in turn, are defined as the

gradients of functions describing the potential energy interactions between atoms of

given types. In this section, the functional forms governding the relevant interactions

are presented.

2.2.1     Potential Functions for Polystyrene

   As is the case with many model materials, the potential energy arising from the

interactions of polystyrene atoms is modeled as the sum of bonded and non-bonded

interactions,


   UP S = Ubonded + Unon−bonded = Ustretch + Ubend + Utorsion + Uvdw + Ucoulomb   (2.1)


where the bonded interactions can be further divided into energetic contributions due

to bond stretching, bond bending, and bond torsion, and the non-bonded interactions

can be separated into van der Waals pair interactions and Coulomb electrostatic

interactions. Two models for polystyrene of the form of Equation 2.1 are presented

in the remainder of this section.




                                          10
2.2.1.1   The All Atom Representation

   The all atom (AA) representaion of polystyrene is a molecular dynamics model

that explicitly accounts for all the atoms that constitute the chain molecules. Fig-

ure 2.1(a) shows the AA representation of a polystyrene monomer. For simulations




                 (a) All atom representation    (b) United atom representa-
                                                tion



                          Figure 2.1: Polystyrene monomer




employing the AA representation of polystyrene, a modification of the potential form

developed by Han and Boyd is used.[10] This potential is modified to include explicit

charges of 0.06e on the hydrogen atoms, and corresponding negative charges on the

carbon atoms. Including these small charges should have a negligible impact on the

polymer properties, as noted by Kaminski.[11] In addition, the bond length of all C-H

bonds is fixed using the SHAKE algorithm[12] implemented in LAMMPS.[7] Again,

this modification is assumed to have a minimal impact on the dynamics of the chain

molecules as far as properties such as Tg are concerned.[13] Finally, the actual values


                                               11
obtained by spectroscopy are used in place of the reduced values used by Han and

Boyd. The potential function is given by

      UP S−AA = kB (l − l0 )2 + kθ (θjik − θ0 )2 +              n
                                                               kφ [1 + dn cos (nφijkl )]
                                                     i=2,3,6
                                                                                           (2.2)
                  qi qj           σ   12     σ   6
               +          + 4ǫ             +
                 4πǫ0 rij         r          r

where rij is the distance between the i-th and j-th atoms, θjik is the angle formed

by the atom triplet (j, i, k), and φijkl is the dihedral angle of the atom quadruplet

(i, j, k, l). The values of the coefficients are given in Table 2.1.

   Due to the computational load required to simulate a system based on this model,

all-atom polystyrene is used only as a tool for comparison. In this work, a reduced

order model is employed in all simulations from which data is collected. The behavior

of the reduced order model is found to approximate the behavior of the corresponding

real systems quite well.

2.2.1.2   The United Atom Representation

   A united atom model of the structure that invokes coarse-graining to reduce the

degrees of freedom is shown in Figure 2.1(b). In the united-atom representation, each

polystyrene monomer is represented by eight united atoms, C1 to C8 , of five different

types. The two outer backbone units for each PS chain, CH3 , and the inner backbone

unit, CH2 , are represented by C1 . The other inner backbone unit CH is represented

by C2 in the united atom model. The C3 unit corresponds to the aromatic carbon

atom without the hydrogen (Caro ) and the units C4 to C8 represent the aromatic

carbon atom with hydrogen (CHaro ) in the phenyl ring.

   The effect of coarse-graining an alkane chain with a united atom (UA) represen-

tation is studied in [14, 15] and it is shown that the C-H bond atoms can be safely

                                            12
coarse-grained into united atoms for the study of the properties of pure polystyrene.

The suppression of the C-H bond vibration in the UA scheme due to united-atom rep-

resentation allows for the time-step of a simulation to be increased by a factor of two

(see Appendix B). Moreover, there is a reduction by half in the number of degrees




Table 2.1: Parameter values (energies in kcal/mol) in the potential energy functions
for the all atom model of polystyrene. The C-H LJ parameter is calculated via the
Lorentz-Berthelot mixing rule.[1]


                             Bond stretch parameters
                  Bond type                    kB            r0 [˚]
                                                                 A
                  C-C                        316.770           1.54
                  C-H                        316.770           1.09
                  C-Caro                     525.561           1.50
                  Caro -Caro                 316.770           1.39

                                 Bond angle parameters
                  Bond type                      kθ         θ0 [◦ ]
                  C-C-C                         57.573       111.60
                  C-C-H, C-C-Caro ,             57.573       109.47
                  Caro -C-H, C-C-Caro
                  Caro -Caro -Caro              71.906       120.00

                                Bond torsion parameters
                                     n
                  Bond type         kφ           dn            n
                  -C-C-             1.6006       1             3
                  -C-Caro -         0.2461       1             6
                  -Caro -Caro -    12.9001      −1             2

                            LJ pair potential parameters
                  Interaction type                ε          σ [˚]
                                                                A
                  C-C                            0.09484      3.450
                  H-H                            0.00979      3.000
                  C-H                            0.03048      3.225



                                          13
of freedom in the system, so that effectively four times the amount of polystyrene

can be simulated with the same computational resources. Given these advantages,

the UA model is preferred where measurements of the quantities of interest are not

affected by this reduction in the order of the model. The potential function for UA

polystyrene is given by
                                                 3
                                       2    1          j
                 UP S−U A   =kθ (θ − θ0 ) +           kφ cos 1 + (−1)j+1 jφ
                                            2   j=1                               (2.3)
                               1                          σ   12     σ   6
                              + kψ (ψ − ψ0 ) + 4ǫ                  +
                               2                          r          r

where the parameters found in Table 2.2 are given in the papers by Wick et. al.[16] and

Harmandaris et. al.[17] Actual values from spectroscopy are used for the parameters.

   One of the artifacts of coarse graining is the possibility of loss in this configu-

rational integrity due to removal of the hydrogen atoms. For example, in the AA

polystyrene model, the planarity of the phenyl ring is maintained via dihedral barri-

ers for the Caro -Caro -Caro -H angles. In the UA representation, a combination of OPLS

dihedral and harmonic improper bending potentials are used to maintain planarity of

the phenyl ring. Additional harmonic improper potentials are introduced to maintain

the conformation of the united atoms along the chain backbone.

2.2.2    Potential Functions for Carbon Dioxide

   Carbon dioxide is modeled as a linear three atom molecule with two oxygen atoms

sharing double bonds with a carbon atom in the middle, as shown in Figure 2.2. The

elementary physical model (EPM2) potential developed by Harris and Yung [18] is

used to model CO2 . The bonded potential consists of contibutions from harmonic

bond stretching and bond bending. The non-bonded van der Waal and coulombic


                                            14
Table 2.2: Parameter values (energies in kcal/mol) in the potential energy functions
for PS. The mixed LJ parameters are calculated via the Lorentz-Berthelot mixing
rule.[1]


                              Bond stretch parameters
                    Bond type                     kB        r0 [˚]
                                                                A
                    CHx -CH                     316.92       1.54
                    CH-Caro                     316.92       1.51
                    C(H)aro -CHaro              525.81       1.40

                               Bond angle parameters
                    Bond type                     kθ        θ0 [◦ ]
                    CH3 -CH-CH2                  62.14       112
                    CH-CH2 -CH                   62.14       114
                    C(H)aro -C(H)aro -CH(aro)   119.50       120

                             Bond torsion parameters
                                          1        2           3
                    Bond type            kφ      kφ           kφ
                    CHx -CH-CH2 -CH     1.41    -0.27        3.14

                       Bond improper dihedral parameters
                    Bond type                     kψ   ψ0 [◦ ]
                    C(H)aro -CHaro -CHaro -CHaro 20.00   0
                    CH-CHx -CHy -Caro            40.01 35.26

                             LJ pair potential parameters
                    Interaction type                  ε     σ [˚]
                                                               A
                    CH3 -CH3                       0.1950   3.750
                    CH2 -CH2                       0.0915   3.950
                    CH-CH                          0.0199   4.650
                    Caro -Caro                     0.1003   3.695
                    CHaro -CHaro                   0.0596   3.700




interactions between atoms belonging to different CO2 molecules are modeled by LJ-

type and charge-charge interactions, respectively. The full potential is given as
                          1                1                 qi qj
                  UCO2 = kB (r − r0 )2 + kθ (θ − θ0 )2 +
                          2                2               4πε0 rij
                                          12           6
                                     σij 15      σij
                            + 4ǫij          −
                                     rij         rij
                             Figure 2.2: CO2 molecule




where the values of the parameters are given in Table 2.3.




Table 2.3: Parameter values (energies in kcal/mol) in the potential energy functions
for CO2 . The mixed LJ parameters are calculated via the Lorentz-Berthelot mixing
rule.


                            Bond stretch parameters
                       Bond type            kB      r0 [˚]
                                                        A
                       C=O                557.45    1.149

                             Bond angle parameters
                       Bond type           kθ      θ0 [◦ ]
                       O=C=O             147.71     180

                        Coulombic interaction parameters
                       Atom type           q [e]
                       O                  −0.3256
                       C                  +0.6512

                           LJ pair potential parameters
                       Interaction type       ε      σ [˚]
                                                        A
                       O-O                 0.1937    3.033
                       C-C                 0.1228    2.757
                       C-O                 0.1542    2.892



                                         16
   Validation of the parameters for bulk CO2 simulated with EPM2 can be found in

literature. [19] As an additional validation in this study, the diffusion coefficient is

obtained from MD simulations using these parameters. The corresponding value is

2.17 cm2 /s, which is in good agreement with the value of 2.023 cm2 /s, obtained from

NMR spin-lattice relaxation time measurements.[20]

2.2.3    Potential Functions for Silica

   Silica is modeled in its amorphous form with a stoichiometric ratio of two oxygen

atoms for every silicon atom. The stucture of amorphous silica is locally tetrahedral,

with a unit forming as shown in Figure 2.3. The bulk material is composed of twofold-




                         Figure 2.3: Representation of silica.




coordinated oxygen atoms and fourfold-coordinated silicon atoms arranged into a

network of tetrahedra that share corners.[21] The short-range order of the material is

very high, but the repeating structure breaks down rapidly, disrupting the regularity

of the pair-correlation funciton, and making the probing of the detailed structure with

diffraction experiments difficult.


                                          17
   The most popular potential for the MD simulation of amorphous silica is due to

van Beest, Kramer, and van Santen.[22] While there exist more complex potentials

involving multibody interactions,[23, 24, 25] the BKS potential is employed for its

computational efficiency and ability to reproduce the structual features of silica to

the required degree of accuracy. A modified form of this function due to Singer is used

in this work.[26] As shown in Figure 2.4, the modification alleviates the unphysical


                                          250
                                                                                                       BKS
                                                                                                       Patched
                                          200



                                          150
                 En er gy [kca l/ mo l]




                                          100



                                           50



                                            0



                                          −50
                                            0.5   1     1.5   2    2.5            3   3.5        4     4.5       5
                                                                         r [˚ ]
                                                                            A




Figure 2.4: A plot of the Si-O pair potential energy for the modified BKS potential.




negative well that occurs in the orginal equation by ’patching’ the potential for small

values of r:
                                                                          r                 C6
                                                       qi qj      A exp − ρ −               r6
                                                                                                     r < r0
                                          USiO2 =            +                                                       (2.4)
                                                      4πǫ0 r      B+       C12
                                                                                                     r ≥ r0
                                                                           r 12

The values of the coefficients are given in Table 2.4.

   In the BKS potential, there are no bonded interactions. The potential is simply a

combination of electrostatic interaction with the exponential repulsion and 6th power

                                                                   18
Table 2.4: Parameter values (energies in kcal/mol) in the potential energy functions
for SiO2 . Note that there is only a coulombic part to the Si-Si interaction.

                                Coulomb Parameters
                                Atom
                                 type       q [e]
                                  O         -1.2
                                  Si        +2.4

                             BKS potential parameters
     Atom
     Pair         A              ρ [˚]
                                    A         C6 [˚6 ]
                                                  A          B      C12 [˚12 ]
                                                                         A
     O-O        32025.8680         0.3623     4035.5961     91.2495   673.8907
     O-Si      415176.5281         0.2052     3079.4619    107.6398 11818.2927




attractive terms of the Born-Huggins-Mayer (BHM) potential [27] to represent the

van der Waal interaction. This potential has been found to reproduce many important

structural features of silica across a range of temperatures. [28] To implement the

BKS function in LAMMPS, it was necessary to derive the force due to the potential,

which is

                                               A        r     6C6
                       ∂         qi qj         ρ
                                                 exp   −ρ −    r7
                                                                    r < r0
            FSiO2   = − USiO2 =          +                                       (2.5)
                       ∂r       4πǫ0 r 2       12C12
                                                                    r ≥ r0
                                                r 13

   The surface of a slab of amorphous silica requires special attention. When exposed

to moist air, the slab surface undergoes a hydroxlation process, by which OH groups

are formed in regions containing specific heterogeneities. The modeling of this process

is discussed in some detail in Section 5.2.1. The introduction of hydrogen and new

types of binding neccessitate the implementation of a more complicated potential.

In this work, a potential developed by Hassanali and Singer is employed.[9] The



                                         19
functional form is that of a truncated vessal (TV) potential:
                                                  a
    UT V =k (θjik − θ0 )2 θjik (θjik + θ0 − 2π)2 − π a−1 (π − θ0 )3 − DT V
                           a
                                                  2
                     8     8
                    rij + rik                                                             (2.6)
          × exp −
                        ρ8
where rij and rik are the distances between atoms i and j and atom i and k, respec-

tively, and θjik is the angle formed by the atom triplet (j, i, k), given by

                                                   rij · rik
                                 θjik = cos−1                                             (2.7)
                                                    rij rik

where rij and rik are the vectors from atom i to atoms j and k, respectively. The

values of the coefficients are given in Table 2.5. The full potential for the hydroxylated




Table 2.5: Parameter values (energies in kcal/mol) in the TV potential energy func-
tions for the silica slab. Here, OS refers to oxygen atoms in a silanol group.

                              TV potential parameters
      Interaction type      k      θ0 [◦ ]      a                      ˚
                                                                    ρ [A]   DT V    n
      Si-OS -H              3.801   64.284     −0.4237               1.85   28.60   8.0
      H-OS -H             299.787    0.0        0.0                  1.4     0.0    4.0
      OS -H-OS            299.787    0.0        0.0                  1.4     0.0    4.0
      Si-O-Si             691.817    0.0        0.0                  1.6     0.0    4.0




slab combines the energies of the BKS potential and the truncated vessal potential.

   To implement the truncated vessal potential in LAMMPS, it was neccessary to

derive the force contribution from this energy. The following derivation follows the

development in the DLPOLY user manual.[29] As a first step, it is noted that there

is a multiplicative decomposition of the potential


                         UT V = A (θjik ) S1 (rij ) S2 (rik ) S3 (rjk )                   (2.8)

                                              20
The force can then be found by differentiating

                         ∂            ∂
              flα = −        UT V = − α [A (θjik ) S1 (rij ) S2 (rik ) S3 (rjk )]       (2.9)
                        ∂rlα         ∂rl

Applying the product rule and using the definitions
                                       8    8
                                      ri + rik     1       ∂A (θjik )
                         γ = exp −                                                    (2.10a)
                                         ρ8    sin (θijk ) ∂θjik
                                         7
                                        rij      r8 + r8
                        γa = 8A (θjik ) 8 exp − i 8 ik                                (2.10b)
                                        ρ           ρ
                                         7        8      8
                                        rij      ri + rik
                        γc = 8A (θjik ) 8 exp −                                       (2.10c)
                                        ρ           ρ8

it can be shown that the forces on atoms i, j, and k are given by
              α     α                     α            α                α        α
             rij + rik                  rij           rik              rij      rik
  fiα= −γ              − cos (θjik )          +                 − γa       − γc       (2.11a)
             rij rik                    rij 2        rik    2          rij      rik
               α                     α                α
   α         rik                    rij              rij
  fj = γ             − cos (θjik )          + γa                                      (2.11b)
          rij rik                  rij 2             rij
               α
   α         rik                    rα                α
                                                     rik
  fk = γ             − cos (θjik ) ik 2 + γc                                          (2.11c)
          rij rik                  rik               rik

Based on this definition of the forces, the contribution to the virial is


                                 W = − (rij · fj + rik · fk )                          (2.12)


and the contribution to the atomic stress tensor is

                                                      α β
                                    σ αβ = rij fjβ + rik fk
                                            α
                                                                                       (2.13)

2.3     Modeling Component Interactions

   In order to accurately model interface effects, careful attention must be paid to the

interactions between atoms belonging to different component materials. In this sec-

tion, the development of the polymer-gas and polymer-substrate interactions relevant

to this study is presented.

                                              21
2.3.1                               Polystyrene - Carbon Dioxide Interactions

      Interactions between atoms belonging to polystyrene molecules and those belong-

ing to CO2 molecules govern important physical properties of the binary polymer-gas

systems being studied. The non-bonded interactions between these pairs of atoms

are modeled using a LJ potential.[30] For non-bonded interactions between dissimilar

atoms (i = j), the energy well depth (ǫij ) and zero-potential distance (σij ) parame-

ters are typically determined using the Lorentz-Bertholet (LB) mixing rule.[1] In this

method, σij is calculated as the arithmetic mean of σii and σjj , and ǫij is calculated

as the geometric mean of ǫii and ǫjj .

      As shown in Figure 2.5, MD simulations in which the polystyrene-CO2 interactions


                                                                    LB                                                                                LB
                               15                                                                                15
                                                                    CM−HHG                                                                            CM−HHG
                                                                    PALS                                                                              PALS
   Probability Density of fv




                                                                                     Probability Density of fv




                               10                                                                                10




                               5                                                                                  5




                               0                                                                                  0
                                0    0.05   0.1        0.15   0.2        0.25                                      0   0.05   0.1        0.15   0.2        0.25
                                                  fv                                                                                fv

                                            (a) 320K                                                                          (b) 370K



Figure 2.5: Comparison of fv data from simulations of a binary polystyrene-CO2
system with PALS results [2] at ambient pressure



determined via the LB mixing rule do not yield satisfactory agreement with experi-

mental free volume (fv ) results of Jean [2] across the required range of temperatures

for this study.

                                                                                22
   In the simulations, the free volume fraction is measured as the ratio of the accessi-

ble volume to the total volume, where the accessible volume is the difference between

the measured unoccupied volume and the theoretical limit for the unoccupied volume

in a system of randomly packed spheres. The unoccupied volume is measured as the

difference between the total simulation volume and the space filled by atoms, where

each atom is manifested as a sphere with a radius equal to its van der Waals radius,

[8] e.g., those given in Table 2.2. The volume determining algorithm uses cubic voxels

with a resolution of ∼0.01 ˚.
                           A

   In addition, the density measurements for the simulations are approximately 10%

higher than the bulk density value found in experiments. This discrepancy has also

been noted in a paper by Halgren,[31] where the LB mixing rule was found to under-

estimate σij and overestimate ǫij for pair interactions between atoms with significant

differences in the values of these parameters. Thus, it is to be expected that the

density value is overestimated.

   A different mixing rule, previously employed for inert gases,[31] is chosen to repre-

sent the PS-CO2 interatomic non-bonded potential in this study. This rule, denoted

here as ’CM-HHG’, employs the cubic mean for calculations of σij and the harmonic

mean of the harmonic and geometric means for ǫij :

                                                 3   3
                                              σii + σjj
                                        σij =    2   2
                                                                                (2.14a)
                                              σii + σjj
                                            4ǫii ǫjj
                                  ǫij = √        √    2                         (2.14b)
                                          ǫii + ǫjj

The values of the parameters calculated via the CM-HHG rule are found in Table 2.6.




                                           23
   The CM-HHG mixing rule mirrors the procedure used by Sanchez and Lacombe,

[32] wherein the LB mixing rule is ammended with empirical factors ξ and ζ for σij

and ǫij , respectively. This is found to more accurately calculate thermodynamical

properties like heat of mixing, critical solution temperature and chemical potential

using Lattice-Fluid theory for polymer solutions. Analogues of the factors ξ and ζ

can be calculated for the CM-HHG mixing rule:
                           CM
                         σij −HHG           σii σjj
                      ξ=       LB
                                  =2 1− 2           2
                                                        ≥ 1.0                (2.15a)
                             σij           σii + σjj
                                       √ √
                         ǫCM −HHG
                          ij         4 ǫii ǫjj
                      ζ=          = √       √         ≤ 1.0                 (2.15b)
                             ǫLB
                              ij    [ ǫii + ǫjj ]2

Thus, the CM-HHG mixing rule gives higher values of σij and lower values of ǫij than

that obtained by the LB mixing rule, which is consistent with previous work.[32,

33]. The improvement in agreement with experimental measurements of fv is seen in

Figure 2.5.




Table 2.6: Mixed LJ parameters for PS-CO2 interaction calculated via the CM-HHG
mixing rule. Energies are in kcal/mol and distances are in angstroms.

                           LJ pair potential parameters
                           Atom type      ε        σ
                           CH3     O 0.1944 3.4665
                           CH2     O 0.1286 3.6099
                           CH      O 0.0456 4.1674
                           Caro    O 0.1357 3.4285
                           CHaro O 0.0986 3.4319
                           CH3     C 0.1527 3.4016
                           CH2     C 0.1054 3.5592
                           CH      C 0.0404 4.1576
                           Caro    C 0.1107 3.3595
                           CHaro C 0.0828 3.3633


                                        24
2.3.2     Polystyrene - Silica Interactions

   The interaction between atoms in polystyrene molecules and atoms in a silica

substrate has an important impact on the local dynamic behavior of the polymer near

the interface. To accurately model this interaction, a pair potential was developed

based on the ab initio work of Singer.[34] This potential mixes the ab initio determined

LJ paramters of the silica atoms with the polystyrene parameters using the CM-HHG

mixing rule. The values of the parameters are given in Table 2.7.




Table 2.7: LJ parameters for PS-Silica interactions. Energies are in kcal/mol and
distances are in angstroms. The atom types Os and Ox refer to atoms on the surface
and within the silica slab, respectively

                             LJ pair potential parameters
                             Atom type       ε       σ
                             CH3     Os 0.26341 3.04
                             CH3     Si   0.3377 3.568
                             CH3     Ox 0.78375 3.179
                             CH2     Os 0.18044 3.14
                             CH2     Si 0.23132 3.668
                             CH2     Ox 0.53687 3.279
                             CH      Os 0.08415 3.49
                             CH      Si 0.10788 4.018
                             CH      Ox 0.25037 3.629
                             Caro    Os 0.14563 3.015
                             Caro    Si 0.18669 3.541
                             Caro    Ox 0.43329 3.154
                             CHaro Os 0.18892 3.012
                             CHaro Si 0.24219 3.541
                             CHaro Ox 0.5621 3.152




                                          25
2.4    Creation of Polystyrene Initial Configurations

   Before MD simulations can be carried out, simulation volumes must be populated

with atoms of the materials being modeled. It can be a challenging task to do this

in a manner that provides a reasonably stable intial configuration, especially for the

long chain molecule polymer configurations used throughout this work. The method

developed to populate the various simulation volumes with PS is presented in this

section.

   There exist a number of polymer sample preparation methods for molecular dy-

namics modeling, including the Rotational Isomeric State (RIS) model [35, 36, 37]

and the Phantom Chain Growth (PCG) model. [38, 39] These methods are intended

to create polymer chains within a simulation volume with a bias toward expected

equilibrium conformations. In the PCG method, monomers are introduced based on

a random distribution of the backbone dihedral angles. The change in energy ∆U,

resulting from the introduction of a new site is determined by summing all the di-

hedral and non-bonded interactions. The acceptance criterion is a probability p for

accepting the new site

                                 p = min 1, e−β∆U                               (2.16)

where β = 1/ kb T , kb is Boltzmann’s constant and T is the temperature of the system.

If a new monomer segment is not accepted, the chain is shortened by a single segment

and the procedure is repeated again, sampling another dihedral angle.

   To avoid the phenomenon of ’deadlocking’ that can negatively affect the PCG

algorithm during the creation of a large system in a fixed simulation volume, an

augmented PCG scheme is employed in this work. The probability of deadlocking



                                         26
is high for polystyrene, where large phenyl groups are attached to the backbone

carbon atoms, and the presence of boundary atoms adds to this phenomenon. In

the augmented scheme, the acceptance criterion given in Equation 2.16 is relaxed by

calculating ∆U without considering the phenyl side group. If the additional backbone

atom is acceptable to the chain, the corresponding phenyl group is forced onto the

chain at a later stage. This can create an unphyscially high energy locally within the

initial configuration, and must be taken care of prior to equilibration.

   To eliminate the numerical instabilities due to the strong repulsion at distances

rij << σij , a gradual introduction of excluded volume effects must be achieved

in the equilibration procedure. The first method of this introduction is through a

slow increase in the energy of interactions via a regularization based force capping

method.[40, 8] In this method, the non-bonded LJ potential is modified to introduce

force capping,
                                               LJ
                               U LJ (r ) +   dUnb (r)
                  f c−LJ         nb    c        dr
                                                              (r − rc ) 0 ≤ r ≤ rc
                 Unb     (r) =                           rc                          (2.17)
                               U LJ (r)                               r > rc
                                 nb


In Equation 2.17 the r −12 repulsion in the standard LJ potential is replaced by a

much weaker, linearly increasing repulsion for r ≤ rc . From Equation 2.17, it can

be seen that the maximum force, occurring at r = 0, is determined by the value of

rc and increases as rc decreases. Correspondingly, the rc is initially set at 0.9σ and

slowly reduced across a number of pre-equilibration simulation runs until the system

is numerically stable enough for the full potential to be implemented.




                                              27
   Chapter 3: INITIAL INVESTIGATIONS: INTERFACE
      EFFECTS AND GLASS TRANSITION IN FINITE
       POLYSTYRENE-CARBON DIOXIDE SYSTEMS




3.1    Inroduction

   Carbon dioxide is a widely used processing solvent in the polymer industry [41, 42].

The strength of CO2 as a processing solvent lies in its high diffusivity and liquid like

solubility. [43, 44, 45, 46] Adsorbed CO2 is known to depress Tg by more than

50K.[47, 48] This reduction in Tg facilitates processing of polymer at near room-

temperature conditions, which is especially important for nanoscale biomedical appli-

cations. Furthermore, even at relatively high pressures, CO2 has very little impact on

bioactivity when compared to the effects of elevated temperature,[49] which makes it

an ideal agent for polymer-biomolecule systems.

   To take full advantage of the desirable characteristics of the polymer-gas interac-

tion, a detailed understanding of the precise effects of CO2 on polymer chain dynamics

is neccessary. For example, the depression in Tg is more pronounced at the poly-

mer surface, potentially allowing for the modification of the polymer surface while

leaving intact geometric features with sizes much larger than the surface-affected

depths.[50, 49, 4] Positronium lifetime annihilation spectroscopy (PALS) experiments

                                          28
have been used to probe the morphological changes induced in polystyrene, and have

shown that the depression in Tg is accompanied by a significant increase in free vol-

ume. Thermodynamic models have also been used to understand the effects of the

polymer-gas interaction, and when properly calibrated can accurately predict reduc-

tions in Tg for polymers with increasing pressure.[3, 51, 52, 53, 54] However, there

remains a lack of understanding about the mechanisms of the observed phenomenon

that can be surmounted by deeper insight into the atomic-level interactions that

govern the depression in Tg .

   In this study, MD simulations of PS-CO2 systems are performed to probe the

behavior of the polymer in the presence of CO2 at the atomic scale. MD simulations

of binary gas-polymer systems have been performed with the objectives of studying

gas diffusion and solubility in polymers.[10, 55, 56, 57] However, there is a scarcity of

examples in the literature using MD to model the dynamic behavior of the polymer.

                                                          u
One MD model for PS-CO2 systems, developed by Eslami and M¨ ller-Plathe,[57]

predicts Tg for bulk polystyene and analyzes the sorption of CO2 in PS. This model

does not accurately predict the swelling of PS in the presence of absorbed CO2 due

to its limitations in PS-CO2 interaction parameters, ensemble representation and

applied boundary conditions. These limitations are overcome in the present study

through MD modeling of a nanoscale polystyrene-CO2 binary system with a finite

system geometry described in Section 3.3.0.1.

   Results of the MD simulations in this study are validated against PALS experi-

mental observations and the outputs of the thermodynamic models cited above. The

simulation results provide an enhanced understanding of the effect of CO2 on the




                                          29
surface-affected region of the finite PS system. In Section 3.2, a brief survey of exper-

imental methods is presented, including a brief description of the PALS experiment

that provides the framework for the design of the current MD model. The MD sim-

ulation model is explained in Section 3.3. This section includes initial system setup,

potential functions for the materials involved and the pre-simulation equilibration

protocol to prepare the sample for production runs. The results are reported and

discussed in Section 3.4, with concluding remarks in Section 3.5.


3.2    Experimental Methods For Determining Glass Transi-
       tion Temperature

   Experimental methods like dilatometry and differential scanning calorimetry (DSC)

provide good results for bulk polymer systems, but PALS is now established as

one of the most effective experimental methods for probing nanoscale polymeric

materials.[2, 58, 59] The small size of the positronium probe, at ∼ 1.6˚, makes it
                                                                       A

particularly effective in detecting and measuring small free volume cavities. Not only

can PALS measure cavity size and the free volume fraction (fv ), it can also provide

information on the distribution of cavities in the material in the 1-10 ˚ range.[60, 61]
                                                                        A

   The studies of Jean and co-workers and have shown that Tg decreases with in-

creasing pressures of CO2 , with approximately 50K difference in Tg between ambient

and 400 psi CO2 conditions. In the analysis of the PALS results, the drop in Tg is

attributed to the solubilization of CO2 , which acts as a plasticizing agent, creating

a substantial amount of free volume. However, these experiments are not able to

clearly discern the underlying molecular-level physics involved in this Tg depression.




                                          30
The MD model in this study provides the framework needed to explain the experi-

mental observations by facilitating study of the atomic-scale interactions that govern

the Tg effects.


3.3       System Parameters and Simulation Protocol

   In this study, MD simulations are carried out on model PS systems immersed in

CO2 environments corresponding to a range of applied gas pressures. As noted in

the Section 3.1, the application of a constant virial pressure on a periodic simulation

cell is not feasible for the type of systems and process being modeled. An assessment

of methodologies for pressure application for finite systems in the absence of an ap-

plied constant virial pressure is given by Baltazar et. al.[62] In the current work, an

additional challenge is that the pressure-inducing agent (CO2 ) is also a solvent. Con-

sequently, a novel PS-CO2 simulation model is introduced, wherein CO2 acts both as

a plasticizing and pressure-inducing agent.

3.3.0.1    Finite Polystyrene Configuration

   In this work, a finite, cubic simulation geometry is introduced to study polystyrene.

This simulation setup involves a cubic simulation volume with repulsive walls within

which a finite polystyrene sample is contained. To create such a system, a number

of polystyrene chain molecules are grown in a cubic volume via the augmented PCG

algorithm described in the previous section. To contain the growing chains within the

prescribed volume, new monomer units are rejected if their coordinates are within a

certain distance of the edge of the simulation volume. An example of such a system

is shown in Figure 3.1. The walls are modeled with a repulsive LJ-type interaction



                                          31
Figure 3.1: A schematic of the molecular dynamics setup for a finite polystyrene
system.



given by
                                                σwall   12
                                Uwall = ǫwall                                    (3.1)
                                                 r

where the values of ǫwall and σwall are chosen for convenience as 10 kj/mol and 5 ˚.
                                                                                  A

The size of the simulation volume is chosen so that the polymer atoms do not interact

with the walls after equilibration.

   The MD model system is composed of eight polystyrene chains of 320 monomers

each and prescribed numbers of CO2 molecules dictated by the various pressure condi-

tions confined within a finite walled system to maintain CO2 pressure. The simulation

volume for the 8x320 system is cubic with a side length of 120 ˚. The polystyrene
                                                               A

simulation system is equilibrated at an elevated temperature of 600K for 5 ns.




                                          32
3.3.1     Carbon Dioxide Bath

   For a given simulation volume, the number of CO2 molecules is determined by

solving an inverse problem for the van der Waals equation of state for a target pres-

sure. The CO2 molecules are positioned in simulation volumes via an algorithm that

minimizes the total potential energy, which is the sum of that due to the repulsive

part of the LJ pair interactions and the energetic penalty due to the wall repulsion

Uwall given in Equation 3.1. Following the initial positioning, the pure CO2 systems

are brought to equilibrium at 300K during an equilibration phase of 100 ps.

   The number of CO2 molecules in the simulation box is initially estimated using

the van der Waals equation of state for real gases, i.e.:

                                    n2 a
                               p+          (V − nb) = nRT                        (3.2)
                                    V2

where p is the pressure, V is the volume of the simulation box, n is the number of

CO2 molecules, R is the gas constant (8.314JK− 1mole− 1), T is the temperature, and

a and b are empirical parameters for the real gas. Equation 3.2 accounts for non-

bonded interactions between atoms in real gases. The values of the constants a and

b for CO2 used here are 0.3643 Jm3 mole− 1 and 4.25x10−5 m3 mole− 1 respectively.[63]

Independent MD simulations are subsequently carried out to attain a virial CO2

pressure inside the simulation box, expressed as:
                                                           
                               1          1
                          p = nRT +              Fij · rij                     (3.3)
                              V           3
                                                 i,j(i<j)


The second term on the right hand side includes the time averaged (marked by the

overline) work done by the force Fij exerted by atom i on atom j due to the dis-

placement by their relative distance rij . The first part of this equation is the ideal

                                            33
                                           0.07
                                                                         Targeted Pressure
                                                                         vdw Equations
                                           0.06
                                                                         Pure CO2 calculation

                                           0.05




                        Pressure (katms)
                                           0.04

                                           0.03

                                           0.02

                                           0.01

                                             0
                                              0   0.5              1           1.5              2
                                                        No. of CO / Box Volume                 −3
                                                                 2                      x 10




                  Figure 3.2: Calibration of CO2 pressure at 300K.



gas contribution. Figure 3.2 shows the results of the pressure calibration, where the

value calculated by Equation 3.2 is observed to deviate from the simulation values.

Using the form of the van der Waals equation, but with parameter values suggested

by the simulation results, the number of CO2 molecules necessary to achieve a desired

system pressure is calculated. Once the number of molecules required to achieve a

given pressure at 300K is determined, the simulation box is populated.

   The two systems are interfaced and equilibrated for an additional 5 ns at 480K.

Figure 3.3 shows one such equilibrated system with 1528 CO2 molecules, correspond-

ing to a CO2 pressure of 400 psi at 300K.

   The effectiveness of the simulation model is established by comparison with experi-

mental free-volume data from PALS. The 8x320 system size is chosen after conducting

a size-effect analysis with respect to fv . As seen in Figure 3.3.1, the 8x320 system

shows excellent agreement with PALS free volume data whereas a smaller 8x80 system

gives a significant error.




                                                              34
Figure 3.3: A schematic of molecular dynamics setup for finite PS-CO2 system with
repuslive walls at the boundary.



   The target pressures in this study are 0.7, 1.4, 2.8, and 5.5 MPa. The pressures

reported in the results are the CO2 virial pressures developed in the combined PS-

CO2 systems. These differ from the target pressures achieved in the pure CO2 systems



                                                     0.2
                                                                                        PALS Experiment
                                                                                        MD: 8x80 system
                                                                                        MD: 8x320 system
                       Free Volume fraction ( f )




                                                    0.15
                                         v




                                                     0.1



                                                    0.05



                                                       0
                                                           280   300   320    340     360   380     400
                                                                          Temperature (K)




Figure 3.4: Free volume fraction at 400 psi as a function of temperature for different
system sizes. Error bars represent the standard deviation of the measurement.



                                                                           35
mainly due to the volume exclusion of the PS. A typical PS-CO2 molecular system

at 400K with an initial pressure of 400psi is shown in Figure 3.3.

   The number of molecules in the box is kept fixed throughout the simulation pro-

cess. The simulation time is too short to allow for much diffusion of CO2 into or

out of the PS sample. Consequently, the volume fraction of CO2 is maintained in-

dependently of temperature. It is assumed that the deviations in the instantaneous

pressures from the initial values do not have a significant effect on the amount of

CO2 absorbed in the PS ensemble. Indeed, examination of the CO2 density profile

within PS shows little fluctuation across the range of temperatures simulated. Table

3.1 lists the systems considered with different number of CO2 molecules for different


                        Pressure(psi)         14.7 400.0 800.0
                     # of CO2 molecules        50  1528 3676


      Table 3.1: Number of CO2 required to achieve given pressures at 300K .




pressures.

   The validated inter-atomic potentials for pure polystyrene and CO2 discussed

in sections 2.2.1.2 and 2.2.2, repectively, are used for these simulations. Specific

parameters and combinations rules are introduced to accurately model the cross-

constiuent interactions between polystyrene and CO2 , as detailed in Section 2.3.1.

The intermolecular interaction model has a strong effect on the Tg of the PS-CO2

system.

   To calculate Tg , finite system MD models such as that shown in Figure 3.3 and

described in Section 3.3.0.1 are simulated with constant volume and temperature in

                                         36
                                           1
                                                 Density @ 323K
                                                 Fitted Form
                                        0.995


                                         0.99




                        density (gcc)
                                        0.985


                                         0.98


                                        0.975


                                         0.97
                                             0      200       400         600   800   1000
                                                                 time (ps)




Figure 3.5: Evolution of density as the PS system is cooled from 348K to 323K for a
duration of 1 ns (cooling rate of 0.05K/ps). The density converges after around 1 ns
and further simulation does not cause any perceivable change in density.




the Evans NVT ensemble. Polystyrene is modeled in a vacuum (0 psi) and in three

different CO2 pressure conditions, viz. 14.7 psi, 400 psi and 800 psi. In order to

represent the range of behavior from liquid to amorphous solid, simulations follow a

protocol in which the polymer is cooled from 480K to 200K at intervals of 20K for

each of the four pressure conditions. Each step in the cooling process is simulated for

1 ns with a time step of 1 fs. This cooling rate is consistent with previous studies by

the authors.[8] Figure 3.3.1 shows the convergence of density in the MD simulation.

The simulation and experimental cooling rates differ by many orders of magnitude.

However, MD simulations of systems such as those studied here are limited in duration

to tens of nanoseconds at most. It has been shown[8] that a stable density for times

accessible by MD is achieved within the simulation time range. Thus, simulations of

longer duration would incur a higher computational cost without any significant gain.




                                                                  37
                                        1.4
                                                                                        0 psi
                                                                                        14.7 psi
                                        1.2
                                                                                        400 psi

                                         1




                        Density (gcc)
                                        0.8

                                        0.6

                                        0.4

                                        0.2

                                         0
                                          0     10        20       30        40      50            60
                                               dial distance from center of the sample ( ˚ )
                                              Ra                                         A




Figure 3.6: Density variation across thickness at 400K and different pressures of CO2



Consequently, each simulated cooling step is run for 1 ns, where the specific volume

is found to converge to a stable value.


3.4    Results and Discussion

   The MD model is used to study the effect of CO2 on the glass transition temper-

ature Tg of polystyrene across a wide range of subcritical CO2 pressures. This study

is intended to identify and understand the underlying physical mechanisms that fa-

cilitate the reduction in Tg . The MD simulation data is processed to analyze physical

properties like free volume fraction fv , density ρ and mobility at each temperature.

For determining gradients in fv and ρ through the cross-section, the simulation volume

in Figure 3.3 is divided into a series of spherical shells. This process helps delineate

surface and bulk regions, necessary to establish bulk and surface properties.




                                                                38
                                                                                                                                       5
                                                                                                                                            260K
                                                                                                                                            340K
                                                                                                                                            420K
                                                                                                                                       4




                                                                                                           (ri (t) − ri (0)) 2 ( ˚ )
                                                                          density ( finite )




                                                                                                                                 A
                                                                          density ( periodic )
                                                                          f ( finite )
                           1.2                                                v
                                                                          f ( periodic )
                                                                              v

                                                                                                                                       3
   Free volume fraction,

                             1
       Density (gcc)


                           0.8

                                                                                                                                       2
                           0.6




                                                                                                                         ia tm
                           0.4
                                                                                                                                       1

                           0.2


                                                                                                                                       0
                             0
                                 5   10   15    20     25    30     35    40        45           50                                     0    10      20      30      40         50   60
                                            dial distance from center ( ˚ )
                                           Ra                           A                                                                     Distance from center (zones) ( ˚ )
                                                                                                                                                                             A

  (a) Density and free volume fraction for bulk                                                                               (b) MSD(25 ps) at various temperatures
  and finite systems at 400K



Figure 3.7: Gradients in density, free volume fraction and mean squared displacement
for PS-CO2 systems at 400 psi.




3.4.1                                Regions of Material Behavior: Bulk and Surface Prop-
                                     erties

   An important outcome of this study is the identification of clearly demarcated

regions for bulk and surface analysis. Physical properties such as density, free volume,

segmental motion across the thickness are studied to gain insight into the polymer

dynamics. The surface layer is identified from the onset of deviation in ρ and fv from

the bulk behavior which is evident in Figure 3.7(a). The surface region is defined to

begin at the location at which ρ decreases to 90% of the value in the stable core region,

and extends until ρ is less than 5% of the core value. Quantification of the extent

and morphological properties of the surface layer in the presence of CO2 can be used

to determine conditions that preferentially influence the surface properties with little

disturbance to the core region. This is very important for processing devices with

                                                                                                      39
precise nanoscale geometries, allowing greater control over the formation and bonding

processes, which in turn allows the reduction in size of these devices to critical scales.

This is discussed in great detail in Chapter 4.

   To analyze atomic mobility, mean squared displacement (MSD) is calculated for

polystyrene atoms belonging to different layers along the radial direction. The MSD

is a measure of the distance an atom or molecule travels in a given time, given by


                              MSD(t) = |x(t) − x(0)|2                                (3.4)


where x(t) and x(0) correspond to the positions of the atoms at time=t and time=0

respectively. Figure 3.7(b) shows the variation in MSD for polymer atoms situated in

different layers for a binary system with a CO2 pressure of 400 psi. For temperatures

of 260K, 340K and 420K, a gradient in RMS displacement appears along the radial

direction similar to those found in the measurements of density and free volume shown

in Figure 3.7. This signifies that the surface atoms are more mobile than those in the

core region. For example, polymer atoms at the surface of a sample at 340K (below

Tg ) are in fact more mobile than atoms at the core of a sample at 420K (above Tg ).

3.4.2      Effect of CO2 on Polystyrene Tg

   The simulation based Tg of the MD systems is evaluated by analyzing polymer

atoms at the cores of the various samples. The extent of the core region is easily

identifiable from the gradients in ρ and fv shown in Figure 3.7(a). At low radial dis-

tances, the density plot shows fluctuations that are indicative of lack of a statistically

sufficient number of molecules in this region to form a representative volume. The

density measuremnt stabilizes to a value similar to those observed in periodic bulk



                                           40
                                   1.3
                                               0.0 psi
                                               14.7 psi




     Specific Volume (cm 3 g− 1 )
                                               400 psi
                                                                                                      400
                                                                                                                                     Wissinger&Paulaitis
                                                                                                                                     PALS
                                   1.2
                                                                                                                                     MD
                                                                                                      380


                                                                                                      360




                                                                                             Tg (K)
                                   1.1

                                                                                                      340


                                    1                                                                 320
                                         200    240       280   320   360   400   440
                                                          Temperature (K)
                                                                                                      300
                                                                                                            0   200    400       600        800       1000
                                                                                                                      CO pressure (psi)
                                                                                                                         2


         (a) Thermal expansion for MD systems                                                          (b) Comparison of Tg measurements



Figure 3.8: In (a), specific volume is plotted against temperature for finite systems
with different CO2 pressures. The glass transition point is identified by the change in
slope (rate of thermal expansion). Standard deviation in measurements is shown by
error bars. In (b), the determination of Tg from MD is compared with theory[3] and
experiment.[4] The error in the MD measurements is bounded by the error shown in
(a).




systems in the range 10-30 ˚. Beyond 30 ˚, the density rapidly decreases, indicat-
                           A            A

ing surface effects of the molecular system. Similar trends are also seen in the free

volume fraction behavior in Figure 3.7(a). A spherical core region with a 20 ˚ ra-
                                                                             A

dius corresponds to the largest volume for which the studied properties are stable

for all temperature and CO2 pressure conditions. Consequently, this region is used

consistently for analysis of the glass transition behavior of the simulation systems.

   The glass transition temperature of polystyrene in the presence of CO2 is evalu-

ated from the specific volume (ρ−1 ) calculations at different temperatures and CO2

pressure conditions. These values are further corroborated by fv calculations. Figure

3.8(a) shows the variation of ρ−1 at CO2 pressures of 14.7psi, 400psi and 800psi. Two

                                                                                        41
important observations are made from these plots. First, for a given temperature,

the specific volume is found to increase significantly with CO2 pressure in a manner

consistent with experimental observations, a feature which is not found in simulations

of periodic systems at constant virial pressures. Second, the temperature at which

the slope of the specific volume curve is discontinuous decreases with increasing CO2

pressure. This indicates that CO2 plasticizes the PS, causing increased mobility and

resulting in lower Tg . In the Figure 3.8(a), the intersection of the fit lines corresponds

to Tg , where there is change in the dependence of specific volume on temperature.

Figure 3.8(b) compares simulation-evaluated Tg to values predicted by the Wissinger-

Paulaitis model [51] and observed in PALS experiments.[4] Good agreement between

all three sets of data is observed, with a maximum difference of ∼ 10% for all pressure

conditions.


3.5    Conclusions

   An experimentally validated molecular dynamics model is developed to predict

the Tg of a nanoscale polystyrene sample in the presence of CO2 . Pressure control in

the computational system is achieved through the number of CO2 molecules within

the finite simulation volume. Surface and bulk regions are identified via gradients

in density and fractional free volume, and the behaviors of each are captured simul-

taneously and analyzed within a single model. Data extracted from the bulk core

region is compared to the PALS-based experimental data for free volume fraction.[4]

The model is able to accurately predict changes in Tg with changing CO2 pressure.

Observations made in the surface region indicate a preferentially softened surface

layer exhibiting a marked increase in segmental chain mobility in a surface region


                                           42
behaving in a fluid-like manner even below Tg . The model and methods developed

in this study can be easily extended to study the behavior of the polymer in a thin

film configuration, as shown in Chapter 4.




                                        43
      Chapter 4: UNDERSTANDING THE EFFECTS OF
  CARBON DIOXIDE ON INTERFACIAL BONDING OF
                    POLYSTYRENE THIN FILMS




4.1    Introduction

   Polymer bonding at the nanoscale is a critically important part of the processing

and fabrication of biomedical nanodevices for lab-on-a-chip and drug delivery appli-

cations. Conventional techniques such as welding at elevated temperatures or using

chemical adhesives may not be practicable for two reasons:


  1. The required precision of the geometric features in the devices is destroyed;


  2. Biomolecules implanted in the polymer surfaces cannot survive excessive heat

      or harsh chemicals.


Carbon dioxide assisted bonding techniques are being developed to yield a process

that is biologically benign and that preserves nanoscale features. Experimental setups

like the one schematically depicted in Figure 4.1 are being used to study the impact

of macroscopic variables like time, temperature, CO2 pressure, and applied force on

the development of adhesive strength during the CO2 assisted bonding process.[64]

Parallel experimental efforts are being made to understand the nanoscale properties

                                         44
that govern the effectiveness of various processing conditions,[49] but a limited under-

standing of the process at the nanoscale is still a barrier to fabrication. Additionally,

some simulation work has been focused on understanding failure in polymer thin films,

for example in the work of Gersappe and Robbins[65]. In this study, MD simulations

are used to access information about material behavior at the atomic level, thereby

providing quantitative analysis of the fundamental mechanisms of the processes of

adhesive bond formation and failure.

   The CO2 assisted bonding process relies on the plasticization of the polymer to

allow diffusion of chain segments across the bonding interface. A substantial body

of experimental, theoretical, and simulation work [66, 67, 4, 68] has established that

CO2 is able to depress the bulk glass transition temperature or Tb of polystyrene
                                                                 g

(PS) by as much as 65◦ C. This phenomenon effectively plasticizes the polymer at

lower temperatures. The softening is characterized by heightened chain mobility

and increased free volume.[68, 69] It has also been established that there exists a

surface region in PS thin films where molecules undergo the glass transition at an

even lower surface glass transition temperatures.[70, 71] These two effects can be

effectively utilized to bond polymeric nanodevices at near room temperatures in the




Figure 4.1: Schematic representation of an experimental setup for the CO2 assisted
bonding process in PS thin films for fabricating nanoscale devices.


                                           45
presence of CO2 , while maintaining precise geometries and dimensionalities of the

constituents.

   In-situ Bragg diffraction and ex-situ SEM microscopy based tests have shown [50]

that the dissolved CO2 provides the mobility for the polymer chains at the interface

needed for welding. Low temperature fusion of polymeric structure using CO2 has

been successfully carried out for sub-micron structures.[49, 72] However, the CO2

enabled processes have been unable to reliably produce satisfactory bonding in devices

with dimensions less than 100 nm. Optimization of CO2 assisted polymer processing

techniques is therefore a considerable design challenge for the fabrication of biomedical

nanodevices by polymer bonding.


4.2    Adhesive Strength Development at a Symmetric Inter-
       face

   Adhesive strength at a symmetric polymer-polymer interface is known to depend

on the amount of diffusion across the bonding interface.[73] As this diffusion process

takes place, a region develops where atoms from both sides persist. It has been shown

in previous studies that the strength of the adhesion at the interface is dependent on

the width of this region.[74, 75] The rate of diffusion depends on the temperature,

external force and the amount of time elapsed. This study shows that gas sorption

as a result of processing in a CO2 environment also impacts diffusion, increasing

the rate at which polymer thin films can diffuse across the interface. This result is

in accordance with the theoretical findings of Talreja et. al.,[76] who used polymer

density functional theory to posit the widening of a low-density interfacial layer with

increasing CO2 pressure.



                                           46
   The adhesive strength at symmetric PS-PS interfaces has been examined in a

number of experimental studies using various methods including T-peel tests and

lap-shear tests.[77, 78, 79, 80]. At the atomic level, adhesive strength is dependent

on the proximity of atoms from molecules on opposite sides of the interface. Beyond

a threshold amount of atomic separation, the attractive forces due to non-bonded

interactions are insufficient to cause adhesion between opposing atomic surfaces at

the interface.

   This paper presents a novel computational method for evaluating adhesive strength

via MD modeling of deformation and failure. This method, which captures data at

finer scales than is possible with current experimental techniques, allows investigation

of the nanoscale mechanisms at work during the macroscopic processes of bonding

and bond failure. The present MD based study shows that the introduction of CO2

increases the local separation between the film surfaces on either side of the interface.

This phenomenon reduces the adhesive strength of the interface due to the increase in

free volume at the interface that is driven by absorption of CO2 . Results of the MD

simulations conclude that the addition of CO2 molecules into the bonding process

produces a competition between one mechanism that serves to facilitate the bond-

ing process and another that hinders it. This suggests that there exists an optimal

amount of CO2 that when introduced into the polymer system will produce the best

bonding. This optimal amount of CO2 is assessed in this study by examining the

functional relationship between bond strength and the amount of CO2 present during

bonding.




                                          47
4.3     Model Development and Simulation Procedures

4.3.1       Overview of Modeling Scheme

   Molecular dynamics models of binary PS-CO2 systems are simulated for under-

standing the operating atomic-scale mechanisms during CO2 assisted polymer bond-

ing. Each system consists of two PS thin films containing a prescribed number of

molecules in a CO2 environment. The system is subjected to a controlled pressure

environment that is applied in the non-periodic thickness direction by a repulsive wall

ineraction given by Equation 3.1. A representative molecular system for this simula-

tion setup is shown in Figure 4.2. Simulations are carried out in a three-stage process.
         nonperiodic direction




                                                                       CO -limiting walls
                                                   PS-limiting walls




                                 (a) Full System                                                 (b) Surface Representation



Figure 4.2: A molecular simulation system for the thin film bonding process with
CPK representation of CO2 : (a) The PS thin films are represented by dynamic bonds
and are shaded differently for identification, (b) The surfaces of the PS thin films are
represented by an interpolation of the van der Waals surfaces.


                                                                                            48
In the first stage, the initial configurations are created and equilibrated as described

in sections 4.3.2.1 and 3.3.1. After equilibration of each component separately and

subsequent equilibration of the full binary system, the bonding process is simulated.

Lastly, the CO2 is eliminated, and simulations are carried out to test the integrity of

the bonding interfaces under tensile loading. All simulations are carried out with the

LAMMPS software package, with an NVT ensemble implementing a Nose/Hoover

thermostat set to 300K with a relaxation time of 10 fs.

4.3.2     Model Materials and Configuration

   The simulated polymer is monodisperse atactic polystyrene with 200 monomers

per chain molecule. This results in a molecular weight of MW =21 kg/mol. This is

well below the critical entanglement molecular weight for PS, Mce =36 kg/mol,[78] thus

avoiding the effects of chain entanglement that are expected to show up only on much

longer time scales than those investigated here.[81] Two film sizes are simulated in

this study. They contain 10 and 15 chain molecules, respectively denoted as ’10x200’

and ’15x200’. The potential functions governing the interactions of the polystyrene

united atoms are those detailed in Section 2.2.1.2.

4.3.2.1   Thin Film Polystyrene-Carbon Dioxide Configurations

   The PS chains are assembled in the form of freestanding thin films, with periodic

boundaries imposed in the two horizontal directions of the cross-section. Nonperiodic,

finite boundary conditions are enforced in the vertical thickness direction. The PS

molecules are assembled using the augmented PCG scheme described in Section 2.4.

In the nonperiodic direction, an additional penalty potential function Up is imposed




                                          49
to contain the growing chains. This potential energy function is given by

                                               σb   12
                                 Up (d) = ǫb                                     (4.1)
                                               d

where d is the distance to the prescribed simulation boundary, and ǫb =10 kcal/mol and

σb =5˚ are empirically determined Lennard-Jones (LJ) parameters. For a monomer
     A

to be added to a growing chain, the acceptance criterion for the energy ∆U + ∆Up

must be satisfied.

   After chain creation, the films are equilibrated at 500K and then cooled to 300K

over a simulation time of 6 ns. Equilibration takes place in a simulation volume

measuring 8 nm by 8 nm in the two periodic directions of the cross-section. At

the end of the equilibration phase, the height of the film in the 10x200 system is

approximately 5 nm and the height of the 15x200 film is 7 nm.

   For CO2 , the EPM2 model developed by Harris and Yung[18] is used with a

harmonic flexible bond potential, as described in Section 2.2.2. The CO2 simulation

volume is populated using the procedure detailed in Section 3.3.1. The dimensions

of the resulting systems are 8 nm by 8 nm in the two periodic directions of the

cross-section and measure 20 nm in the non-periodic film thickness direction.

   To equilibrate the binary systems prior to the bonding simulation, two identical

simulation volumes, each containing a single freestanding PS thin film are interfaced

after translating each film in the non-periodic direction by ±(h/2 + 0.2nm), h being

the film height. This creates a simulation volume containing two identical freestanding

films separated by a gap of at least 0.4 nm in the nonperiodic direction. The two-film

configuration is equilibrated over 100 ps to create a symmetric PS-PS interface. At

the end of this equilibration, the gap between the films is reduced significantly due

to movement of the film mass centers, but the interface between the two films is

                                         50
still marked by a sharp drop in PS density. This is shown in the density profiles of

Figure 4.3(a). A simulation volume containing the equilibrated CO2 sample is then


                                                 Film 1       Film 2                                    1.25              Film 1       Film 2      Combined


                         1
                                                                                                          1




                                                                                     Density (g/ cm3)
    Density (g/ cm3)




                       0.75
                                                                                                        0.75



                        0.5
                                                                                                         0.5




                       0.25                                                                             0.25




                          0                                                                                0
                         −15   −10      −5             0               5   10   15                        −15   −10      −5             0            5        10   15
                                        Dist an ce from Interf ace (nm )                                                 Dist an ce from Interf ace (nm )


                                     (a) Before bonding                                                               (b) After bonding



Figure 4.3: Polystyrene density along the nonperiodic thickness direction for the
15x200, 800psi thin film system.




superimposed on the two-film volume to create the full PS-CO2 binary system. A

soft potential, given in the LAMMPS literature[7]

                                                                                                  πrij
                                                    E(rij ) = A 1 + cos                                               rij < rc                                     (4.2)
                                                                                                   rc

is initially used between PS and CO2 to eliminate any excessive forces due to the

physically unrealistic overlap that may initially occur between the atoms. The system

trajectories are evolved over a short simulation period of 10 ps in preparation for the

bonding simulations.

   For the remaining simulations, a full LJ potential is implemented between PS and

CO2 atoms. The potential parameters are determined using the CM-HHG mixing rule

described in Section 2.3.1. Parameter values calculated by this mixing rule are given

                                                                                51
in Table 2.6. As detailed in Srivastava et. al.[68], this mixing rule is implemented

for interactions between dissimilar species for its ability to produce experimentally

observed free volume characteristics better than the commonly used Lorentz-Berthelot

rule. [1]


4.4     Investigation of Thin Film Behavior: Model Validation
        and Preliminary Study

   Prior to modeling the actual interface bonding process, a study of the thin film

characteristics is performed to validate parameters in the potential functions and the

overall simulation setup.

4.4.1       Film Density Characteristics

   The thin film model for polystyrene is first validated by comparing density at

the film core to the experimental value of 1.05 g/cm3 reported by Wissinger and

Paulitis.[51] As shown in Figure 4.4, the simulation values for both films of pure PS

are within 1% of the experimental value. Near the film surfaces, however, there exists

a gradient in the density. The extent of this surface region is roughly constant for

the two film heights studied. The addition of CO2 extends this region somewhat, for

reasons discussed below. In addition, CO2 causes an expansion in the PS at the film

core. This expansion is quantified in terms of the change in volume for the pure PS

system shown in Figure 4.5. The agreement with the experimental data of Wissinger

and Paulitis [51] confirm the accuracy of the MD simulation model in predicting

swelling due to the presence of CO2 .




                                         52
                                                              No CO           7 MPa CO
                                                                      2                    2
                                      1.25

                                                15x200
                                        1




                     Density g⊳cm 3
                                      0.75



                                       0.5                                    10x200



                                      0.25



                                        0
                                        −6      −4       −2               0            2       4   6
                                                               Location (nm)




Figure 4.4: Polystyrene density along the nonperiodic thickness direction. A smooth
curve is fitted to highlight the change in profiles upon addition of CO2 .




4.4.2    Chain Dynamics in the Thin Films

   In addition to the variation in surface density characteristics, the efficacy of bond-

ing is determined in large part by the chain dynamics that take place at the bonding

interface. Significant studies have been conducted on the surface dynamics of thin

PS films. It has been argued by Kawaguchi et. al.[82] that there exists a gradient in

mobility across the thickness of the film. Kawana and Jones[83] have demonstrated

that there exists a liquid-like surface layer that persists below the glass transition

temperature. The computational methods used in this study allow the investigation

of surface dynamics at atomic resolution and in the absence of the substrate effects

that are necessarily present in any experimental study. To quantify the deviation in

mobility found near the film surface, a relative measure of MSD is defined


                                             MSDrel = MSDlayer /MSDcore                                (4.3)



                                                                53
                                          9
                                               10x200
                                               15x200
                                          8


                                          7




                      Volume Change (%)
                                          6


                                          5


                                          4


                                          3


                                          2


                                          1


                                          0
                                           0     1      2       3         4        5   6   7
                                                            CO 2 Pressure (MPa )




Figure 4.5: Percentage volume change for two PS thin film under various CO2 pressure
conditions (zero pressure corresponds to absence of CO2 ). Error is of the order of the
marker size.




where MSDlayer corresponds to the mean squared displacement for atoms in a given

layer, and MSDcore is averaged over atoms in all layers from the center of the film,

where there is no layerwise gradient in mobility. As shown in Figure 4.6, the model in

this study is able to quantify the effect of the free surface on mobility. It shows that

an enhancement in the mobility extends ∼1 nm into the film. This effect contributes

to a 2-3 fold increase in the mobility of atoms near a free surface in comparison with

those in the film cores. This increase in mobility at the free surface is not present

at the film-film interface during bonding. However, as discussed below, the addition

of CO2 can similarly increase the mobility at the bonding interface, causing the film

surfaces to behave more like a free surface.




                                                                    54
                                                  2.8
                                                         10x200
                                                  2.6    15x200


                                                  2.4




                      MSD Relative to Film Core
                                                  2.2

                                                   2

                                                  1.8

                                                  1.6

                                                  1.4

                                                  1.2

                                                   1

                                                  0.8
                                                     0   0.5      1      1.5     2     2.5    3        3.5   4   4.5
                                                                      Distance from Film Center (nm)




Figure 4.6: Layer-wise mobility, MSDrel , for polymer atoms in a thin film. Error
bars represent the standard deviation in the measurement.



4.5     Numerical Experiments on the Bonding Process

4.5.1    Simulation of Bonding in PS-CO2 Binary Systems

   As described by Prager and Tirrell,[84] bonding at a symmetric polymer-polymer

interface occurs as the interface ’heals’, eventually becoming indistinguishable from

the rest of the material. During this healing process, the adhesive strength of the

interface gradually grows as the material surfaces on opposite sides of the interface

rearrange and the discontinuities in the material properties that define the interface

are eliminated. In this study, three mechanisms of healing are identified and quanti-

fied. First, the discontinuity in material density at the interface is eliminated as the

film surfaces rearrange to conform more closely to one another. Second, groups of

atoms near the film surfaces diffuse across the boundary until the interface between

molecules of different films begins to resemble the interface between molecules within


                                                                                 55
a film. Third, chain segments from molecules belonging to one film begin to penetrate

into the space occupied by molecules of the opposite film. The numerical experiments

performed in this study are designed to study these mechanisms under various CO2

pressure conditions in terms of their impact on bond integrity.

   Using procedures described above, simulations of the bonding process are carried

out for ten model systems, corresponding to two film heights and five CO2 pressure

conditions. The simulations are conducted at a constant temperature of 300K for a

period of 3.0 ns with a timestep of 3 fs. This simulation time is orders of magnitude

lower than the experimental times, during which Fickian diffusion is thought to be

the main driver of the development of interfacial adhesive strength. Thus, these sim-

ulations attempt only to identify the initial healing of the interface at times shorter

than the Rouse time.[85] Additionally, longer simulations of 33 ns are carried out

for two systems of interest in order to analyze the time-dependence of some of the

important results. During the simulations, statistical data is collected and subse-

quently analyzed using the procedures detailed in the remainder of this section in

order to provide a comprehensive quantification of the atomic behavior, critical to

the development of adhesive strength at the bonding interface.

4.5.2    Density Profiles and Diffusion at the Interface

   At the initial stages of the bonding simulation, the discontinuity in interfacial

density seen in Figure 4.3(a) is the most obvious obstacle to the development of

the adhesive strength. This discontinuity is nevertheless bridged rapidly as seen

in Figure 4.3(b), and all that remains is a small, localized separation between the

surfaces of the films at the interface. To identify this material discontinuity, the



                                          56
local separation between atomic surfaces of the films at the interface is measured by

calculating the perpendicular distance between the van der Waals surfaces of surface

atoms in opposite films over a square grid with resolution 1 ˚. A mean value is
                                                            A

obtained by averaging the separation at each grid point.

   During the bonding simulation, the mean surface separation rapidly decreases and

subsequently stabilizes to roughly its final value at times ∼ 300 ps. This is depicted

in Figure 4.7. In Figure 4.8,




                                                 1.2


                                                  1
                      Surface Separation ( ˚ )
                                           A




                                                 0.8


                                                 0.6


                                                 0.4


                                                 0.2


                                                  0
                                                   0   50   100      150      200   250   300
                                                                  Time (ps)




Figure 4.7: Surface separation at the interface versus time for the 15 chain molecule,
200 psi PS-CO2 system. The curve is typical of all simulation systems, except that
the equilibrium separation S0 depends on CO2 pressure, as seen in Figure 4.8.



    the final separation between the film surfaces increases with increasing CO2 pres-

sure. Thus, according to this measure, the interface becomes more pronounced with

increasing CO2 pressure. This phenomenon indicates that at least one effect intro-

duced by the presence of CO2 serves to degrade adhesive strength. The diffusion of



                                                                  57
                                                  10x200
                                                  15x200

                                         1



                                        0.8




                     Separation ( ˚ )
                                  A
                                        0.6



                                        0.4



                                        0.2



                                         0
                                              0       1       2       3         4       5   6   7
                                                                  CO2 Pressure (MPa )




Figure 4.8: Equilibrium surface separation S0 at the interface following bonding sim-
ulations as a function of CO2 pressure. Error bars represent the standard deviation
in the measurement.



chain segments across the interface between the upper and lower thin films is quan-

tified using the average monomer inter-penetration depth X(t), given by
                                                                   ∞
                                                                  z=0
                                                                      z(t)C(z, t)dz
                                                    X(t) =           ∞                              (4.4)
                                                                    z=0
                                                                        C(z, t)dz

where C(z, t) is the monomer concentration profile and z(t) is the distance from a

diffusing monomer back to the interface, discussed in Wool. [81] For the discrete

atomistic system, this measure is calculated as

                                                                        N
                                                                        i=1 zi (t)mi
                                                           X(t) =                                   (4.5)
                                                                       2 N mi
                                                                            i=1


where N is the number of atoms that have crossed the interface, and mi is the atomic

mass of the i-th atom. The factor of 1/2 is added to account for the symmetric

nature of the films across the interface. As shown in Figure 4.9, the value of the

inter-penetration depth increases with increasing CO2 pressure for both film sizes.

                                                                       58
                                     0.5
                                               10x200 (3 ns)
                                               10x200 (33 ns)
                                    0.45       15x200 (3 ns)




                      Depth ( ˚ )
                                     0.4




                              A
                                    0.35


                                     0.3


                                    0.25
                                           0     1       2        3       4     5   6   7
                                                             CO2 Pressure (MPa)




Figure 4.9: Average monomer inter-penetration as a function of CO2 pressure for two
different timescales and film thicknesses. Error bars represent the standard deviation
in the measurement.



Contrary to the trend seen in the mean separation between the atomic surfaces at

the interface, this measure suggests that the bonded interface becomes stronger with

increasing CO2 pressure. Additionally, the measured values of the inter-penetration

depth are shown to increase appreciably for the longer 33 ns simulations, indicating

that a stronger interface should be observed for the resulting systems. Thus, it is seen

that diffusion at the interface is measurably affected upon introduction of CO2 into

the system. However, the overall impact on the development of adhesive strength is

not clearly elucidated by these simple measurements.

4.5.3     Free Volume Characteristics and Chain Bridging at the
          Interface

   To further investigate the atomic structure of film surfaces at the bonded interface,

two additional measures are introduced. First, the free volume at the surface is

characterized in terms of the number and average size of cavities at the interface. A

                                                                 59
cavity is assumed to be a contiguous set of grid points where the separation between

film surfaces is greater than 1 ˚. The cavity volume is calculated as as discrete sum,
                               A

expressed as:
                                                                     N
                                                       VC = A              hi                        (4.6)
                                                                     i=1

where A is the area of a grid square, N is the number of grid squares occupied by

the cavity, and hi is the separation between the film surfaces at the grid square.

As shown in Figures 4.10 and 4.11 respectively, the number of cavities decreases,


                                          40
                                                                                10x200 (3 ns)
                                                                                10x200 (33 ns)
                                          35                                    15x200 (3 ns)
                     Number of Cavities




                                          30

                                          25

                                          20

                                          15

                                          10
                                               0   1   2        3       4     5       6          7
                                                           CO2 Pressure (MPa)




Figure 4.10: Average number of cavities at the interface as a function of CO2 pressure
for two different timescales and film thicknesses. Error bars represent the standard
deviation in the measurement (for the 33 ns simulations, the error is smaller than the
marker size). Lines are a guide to the eye.



while their mean sizes increase with increasing CO2 pressure. The combined effect of

these two trends results in an increasing free volume with increasing CO2 pressure.

This is the main contributor to the mean surface separation trend, detailed above.

With this information, it is possible to reconcile the opposing trends in diffusion


                                                                60
                                             1.4
                                                       10x200 (3 ns)
                                             1.2       10x200 (33 ns)
                                                       15x200 (3 ns)




                     Cavity Volume (nm 3 )
                                              1

                                             0.8

                                             0.6

                                             0.4

                                             0.2

                                              0
                                                   0     1       2        3       4     5   6   7
                                                                     CO2 Pressure (MPa)




Figure 4.11: Average volume of cavities at the interface as a function of CO2 pressure
for two different timescales and film thicknesses. Error bars represent the standard
deviation in the measurement. Lines are a guide to the eye.




observed above. On the one hand, CO2 plasticizes the polymer, enhancing diffusion

across the interface. At the same time, pockets of adsorbed CO2 on the film surfaces

merge and persist, inhibiting close surface contact and increasing the mean surface

separation. It is also noted that changes in the cavity number and volume values are

not statistically significant for the extended 33 ns simulations, indicating that the

free volume characteristics do not contribute significantly to any changes observed in

the systems resulting from the longer simulation time.

   Next, the number of chain segments that bridge the interface is considered. The

existence of a bridging segment is identified by a contiguous set of grid points, where

the penetration of one film into the other is greater than the van der Waals radius of

the PS united atoms. Figure 4.12 shows a remarkable trend in the number of bridges.

It increases with the CO2 pressure reaching a maximum at ∼ 2 MPa and subsequently

decreases with further pressure elevation. The attainment of the plotted values occurs

                                                                          61
                                                   25                                  10x200 (3 ns)
                                                                                       10x200 (33 ns)




                     Number of Bridging Segments
                                                                                       15x200 (3 ns)


                                                   20




                                                   15




                                                   10
                                                        0   1   2        3       4     5     6          7
                                                                    CO2 Pressure (MPa)




Figure 4.12: Number of chain segments bridging the interface as a function of CO2
pressure for two different timescales and film thicknesses. Error bars represent the
standard deviation in the measurement. Lines are a guide to the eye.



in the first few nanoseconds of the simulations, and the values do not change in a

statistically significant way beyond 3 ns for the two extended simulations. Thus,

the development of chain bridges does not contribute significantly to any changes

observed in the systems resulting from the longer 33 ns simulation times. This is the

first direct measurement of a quantity directly related to strength development that

has a maximum at a finite CO2 pressure. The goal of the remainder of this section is

to identify a dynamic mechanism that can explain why this maximum is observed.

4.5.4    Dynamic Behavior of Chains at the Interface

   A study of the dynamic behavior of atoms at the bonding interface sheds light

on the trend in the development of bridging chain segments. The mean squared

displacement for atoms at the interface is measured as given in Equation 3.4. The

data is gathered over a period of 200 ps during the bonding process simulations. In


                                                                         62
a volume with a height of 1 nm centered at the interface, the mobility of atoms is

found to increase with CO2 pressure up to a maximum at ∼2 MPa. Increasing the

CO2 pressure beyond this begins to limit mobility, as shown in Figure 4.13.


                                  0.04                                              10x200
                                                                                    15x200




                                  0.03
                     MSD (nm2 )




                                  0.02




                                  0.01




                                    0
                                         0   1   2        3        4        5   6            7
                                                     CO 2 Pressure (MPa )




Figure 4.13: Mean-squared displacement over 200ps as a function of CO2 pressure.
Error bars represent the standard deviation in the measurement.



    Thus, for a particular range of CO2 pressures, the mobility at the interfacial

surfaces can be enhanced in a manner similar to that observed at the free surfaces

studied. This trend is similar to that found in the number of bridges being developed.

This infers that an increase in mobility promotes chain bridging.

4.5.5    Simulation of Debonding

   After simulating the bonding process, a debonding simulation is performed to test

the adhesive strength of the PS interfacial bond. At the outset of the debonding sim-

ulation, interaction between PS and CO2 atoms is turned off. This enables modeling

of the bonded PS thin films in the absence of CO2 . A layer of atoms 1 nm thick is

                                                         63
frozen at the edges of the films away from the interface in the non-periodic thickness

direction. These atoms are also displaced with a uniform velocity so that the effective

engineering strain rate is e=2e8 s−1 . The bonded films are strained to failure over a
                           ˙

period of 0.5 ns, with failure occurring at ǫf ≈ 6%. This is shown in Figure 4.14. In


                           0.0 MPa     0.8 MPa     1.7 MPa     3.5 MPa      6.8 MPa                                0.0 MPa     1.1 MPa     1.8 MPa     3.3 MPa      6.9 MPa
                   300                                                                                      300



                   250                                                                                      250



                   200                                                                                      200
    Stress (MPa)




                                                                                             Stress (MPa)
                   150                                                                                      150



                   100                                                                                      100



                   50                                                                                       50



                     0                                                                                        0
                      0     0.01       0.02       0.03       0.04        0.05         0.06                     0    0.01       0.02       0.03       0.04        0.05         0.06
                                                 Strain                                                                                  Strain

                                     (a) 10x200 System                                                                       (b) 15x200 System



Figure 4.14: Stress-strain curves for all debonding simulations. Carbon dioxide pres-
sures are given in the legend.




every system simulated, the films separate completely, with failure occurring at the

original interface.

4.5.6                     Interface Integrity: A Strain Energy Density Based
                          Metric

   The importance of the chain bridging and mobility data is demonstrated by the

measurement of integrity of the bonded interface. To measure the bond integrity

in debonding simulations, the stress is calculated for a collection of atoms initially

contained in a volume including the entire cross-section perpendicular to the thickness


                                                                                         64
direction and spanning 2 nm around the interface in the thickness direction. The

stress tensor is calculated from the sum of the local atomic virial stresses as [86, 87]

                       1         1                     ∂Φ rij
                  σ=                           rij ⊗                    − mi vi ⊗ vi   (4.7)
                       V   i
                                 2   j=i
                                                       ∂rij rij

Here, rij is the vector from the position of the i-th atom to the position of the j-th

atom. The velocity vi is the total velocity vector of the i-th atom less the contribution

due to the affine transformation of the deformation. From this measure of stress, the

strain energy is calculated as
                                          ǫf                   ǫf
                               W =             σ : dǫ ≈             σ33 dǫ33           (4.8)
                                      0                    0


Here, ǫ is the Cauchy strain tensor, and ǫf denotes the strain developed at bond

failure, where the maximum normal stress σ33 is achieved. The strain ǫ33 is defined as

the ratio ∆L/L0 of the change of length in the nonperiodic direction to the original

length in that direction. It is assumed that the in-plane stresses and strains are

significantly smaller than those (σ33 , ǫ33 ) in the thickness direction. This measure

of strain energy is used as a metric to determine the relative bond integrity of each

molecular system. Figure 4.15 demonstrates that the interfacial integrity is enhanced

by CO2 pressure conditions up to a maximum at ∼2 MPa. Further addition of CO2

beyond this begins to degrade the integrity of the interface.

   The strain energy density measured for the system consisting of two 10x200 PS

films under 1.7 MPa CO2 is 7.67 MJ/m3 , which is greater than 50% of the value of

14.56 MJ/m3 obtained for a virgin 20x200 PS film subjected to a similar simulation

of failure, indicating a significant development of strength. A sensitivity study is

performed to test the dependence of the strain energy values on time across the range

of MD-accessible simulation times. This study is performed by simulating interface

                                                  65
                                                     8                                         10x200
                                                                                               15x200



                                                    7.5




                     Strain Energy Density ( MJ )
                                             m3
                                                     7




                                                    6.5




                                                     6




                                                    5.5
                                                      0   1   2       3       4        5   6            7
                                                                  CO2 Pressure (Mpa)




Figure 4.15: Strain energy density developed up to failure as a function of CO2
pressure.




healing between two 10x200 films for 0.5, 1.0, 2.0, 3.0, 6.0, and 33.0 ns under no CO2

and 1.7 MPa CO2 pressure conditions and subsequently testing the faliure behavior

as described above. As shown in Figure 4.16 that there is an initial phase of limited

strength development for times less than 1 ns, followed by a rapid development of

strength in the 1-3 ns range, and a final regime of slow strength development beyond 3

ns, continuing to 33 ns. Presumably, this longer-time development continues beyond

33 ns up to experimental time scales.


4.6    Conclusions

   This molecular dynamics based study is intended to investigate the effect of CO2

on the interface of polystyrene (PS) thin films. This property is important in the

fabrication of nano-scale film structures by bonding thin films. Thin film bonding

simulations are carried out on 5 nm and 7 nm thick PS thin films at 300K. This

                                                                      66
                                                    9
                                                           1.7 MPa CO
                                                                    2
                                                           No CO2

                                                    8




                     Strain Energy Density ( GJ )
                                             m3
                                                    7




                                                    6




                                                    5




                                                    4 −1                     0                       1    2
                                                    10                  10                      10       10
                                                                             Simulation Time (ns)




Figure 4.16: Strain energy density developed up to failure as a function of simulation
time for the 10x200 polymer system under two CO2 pressure conditions.




temperature is below the bulk glass transition temperature Tb , but above the surface
                                                            g

glass transition temperature Ts .[71] This temperature is sufficiently high to allow
                              g

for limited diffusion of PS chain segments across symmetric interfaces between the

thin films. The introduction of increasing pressures of CO2 is found to affect certain

important properties of the polymer at the bonding interface, producing a strong

dependence between bond integrity and CO2 pressure.

   Part of this dependence is brought about by two competing trends. First, the

average monomer inter-penetration depth is found to increase with increasing CO2

pressure. This indicates that the introduction of CO2 facilitates the diffusion of chain

segments across the interface, with increasing amounts of diffusion for CO2 pressures

up to ∼ 7 MPa. Second, the introduction of CO2 into the polymer systems is found to

impact the free volume characteristics of the bonding interface. At elevated pressures

of CO2 , free volume cavities at the interface are fewer and larger, and account for


                                                                                   67
a net gain in free volume at the interface. This increase in free volume leads to a

greater average separation between the film surfaces for the range of CO2 pressures

investigated, up to ∼ 7 MPa.

   The competing mechanisms point to an optimal amount of CO2 required for a

maximum enhancement of adhesive strength at the bonded interface. Much of the

strength that is developed is due to the van der Waals attraction between atoms on

opposite sides of the thin film interface. Processing in a CO2 environment facilitates

diffusion of atoms across the interface. It increases the surface area of the films

at the interface and therefore the number of atoms from opposing films that are in

close proximity. This phenomenon enhances the bond strength. At the same time,

CO2 enhances the free volume in the polymer. Separation between the film surfaces

increases, limiting the strength of the attractive forces that atoms at the film surfaces

exert across the interface.

   A second source of adhesive strength in the bonded interface is the existence of

chemically bonded atoms from molecules in one film, bridging the interfacial gap and

penetrating into the opposing film. It is shown that this process is enabled by the

enhanced mobility that is maximized under CO2 pressure conditions in the neigh-

borhood of 2 MPa. As a bridging segment from one film penetrating into the other,

the adhesive strength of this bridge approaches the strength of the covalent bonds

in the segment. Thus, the energy required for debonding at a bridge is potentially

much higher than that required to separate the films when the surfaces are held to-

gether with the weaker van der Waals attraction between the film surfaces. In fact,

increasing the number of bridging segments is demonstrated to be a mechanism for

the development of significant adhesive strength, as previously theorized.[88, 89]


                                          68
   This study finds that the strongest bond is produced for CO2 pressure conditions

in the neighborhood of 2 MPa in the model systems studied. Since the sorption of

CO2 has a negligible affect on the tensile strength[90] and elastic modulus[91] of PS,

it is asserted that the observed increase in bond strength is due to enhancement of the

interfacial adhesion. While the exact values of the temperature and CO2 pressure are

almost certainly material- and size-dependent, two main ideas are widely applicable.

First, for nanoscale polymer structures, the difference between bulk and surface glass

transition temperatures can be exploited to produce viable bonds at a processing

temperature Ts <Tp <Tb . Second, the enhancement of adhesive strength at the
             g       g

bonding interface via CO2 is directly related to the enhancement of chain mobility at

the interface, which is maximum at a finite CO2 pressure.

   Lastly, it is noted that the MD timescale prevents the simulation of the complete

bonding process. To circumvent this limitation, a careful examination is made of

the time-dependence of the quantities investigated to explain the observed trends in

adhesive strength with respect to CO2 pressure. As a result of this analysis, it is

expected that the conclusions made based on the results of this study should remain

valid even for more realistic timescales.




                                            69
      Chapter 5: QUANTIFICATION AND ANALYSIS OF
EFFECTS OF THE SILICA-POLYSTYRENE INTERFACE
      ON THE GLASS TRANSITION IN A SUPPORTED
                     POLYSTYRENE THIN FILM




5.1    Introduction

   Recent advances in techniques for the fabrication of thin polymer films have en-

abled the production of ultrathin films with thicknesses on the order of 10 nm. These

developments are promising for their potential applications in nanoscale biomedical

devices. However, the structural and dynamic properties of these films is signifi-

cantly divergent from the bulk polymer. It has been noted in many experimental and

computational studies that the effects of polymer-substrate interaction and nanocon-

finement contribute greatly to the alteration of certain important features of the

polymer, such as Tg . The Tg of the polymer is of great interest in its own right,

and is also important for its connection to the elastic modulus as noted by Tor-

res et. al.[92] Additionally, an equivalence has been drawn between the behavior of

polymer thin films and polymer nanocomposites,[93] in that both exhibit similar de-

viations in behavior near the polymer-substrate or polymer-particle interfaces.[94] It

has been suggested that the properties of polymer nanocomposites can be fine-tuned

                                         70
by incorporating nanoparticles which have specific interfacial interactions with the

surrounding polymer matrix.[95] Thus, it is seen as an important problem to analyze

and quantify the various deviations in the behavior of the polymer thin films near

such interfaces.[96, 97]

   In this study, MD modeling of a polystyrene thin film supported on a silica sub-

strate is carried out in order to measure the deviation in the local Tg of the polymer

near the substrate. The details of the simulations are layed out in the next section,

followed by a discussion of the results in Section 5.3.


5.2     Simulation Details

   The methods for the simulation of glass transition dynamics are described in this

section. First, initial configurations of silica and various polystyrene geometries are

created and equilibrated. For the binary system, the polymer and substrate systems

are then combined and re-equilibrated. Finally, all the systems are run through the

protocol described in Section 5.2.4.

5.2.1     Generation of Initial Configuration: Silica Slab

   To create a silica slab, a two step process is employed, wherein a periodic simula-

tion volume is populated with a stoichiometric amount of silicon and oxygen atoms

and equilibrated at high temperature. Subsequently, the silica sample is cooled, and

then it is cleaved and equilibrated further at low temperature.

   The population of the initial volume is carried out via an algorithm that places

atoms one at a time, randomly selecting atomic positions in an attempt to minimize

the potential energy of the system, as calculated by Equation 2.4. Following a mod-

ification of the procedures developed by Woodcock et. al. [98] and Soules [99], the

                                          71
initial random atomic configuration is first simulated in a periodic cell with a temper-

ature of 10000K for a period of 100 ps. Subsequently, the ’bulk’ sample is quenched

to 300K over the course of 100 ps. To produce the slab configuration, the bulk sample

is cleaved by removing the periodicity in the z-direction and equilibrating at 300K

for 50 ps.

   The surface of the slab is modified by creating hydroxyl groups at specific defect

sites to create silanol groups. Types of defect include threefold coordinated silicon

sites, two-membered (2M) rings, and nonbridging oxygens (NBO), as shown in Fig-

ure 5.1, and a number of studies have focused on the formation of silanols at these




 (a) Threefold-coordinated Sili-   (b) Two-membered ring    (c) Nonbridging oxygen
 con



Figure 5.1: Pictorial catalogue of SiO2 surface defects. In each figure, the darker
atom(s) constitute the defect. The fictitious ’bonds’ are drawn to highlight the coor-
dination of each atom.




sites, for example Masini nad Bernasconi[100], Hassanali and Singer,[9] and Wilson,

Walsh and co-workers.[101, 102] Following Hassanali and Singer[9], this surface mod-

ification in this work places geminal silanols at NBO sites and vicinal silanols at the

sites of the 2M rings, as shown in Figure 5.2. For the slab prepared in this study,

                                            72
                                                           H           H


                       O      O        O
                   H                       H          O        O           O




                                  Si
                                                      Si                   Si
                                                                   O



               (a) Geminal silanol created at a (b) Vicinal silanol created at
               non-bonded oxygen site           the site of a 2M ring



Figure 5.2: Representation of the creation of silanol groups on the silica slab surface.
The dashed circle represents the original oxygen placement.




this hydroxylation procedure yielded a final silanol density of ∼ 4.66 nm−2 , which is

in line with the experimental results of Zhuravlev[103] and those cited by Iler.[104]

Of these silanols, 26% are vicinal silanols and 74% are geminal silanols.

5.2.2     Generation of Initial Configuration: Polystyrene Thin
          Film

   Using the augmented PCG algorithm described in Section 2.4, polystyrene sam-

ples are produced consisting of 10 molecules of 789 monomers each, resulting in a

63130 united atom simulation system representing polystyrene of molecular weight

MW =82 kg/mol. This atomic configuration is produced in two geometries: 1) Bulk

configuration with period boundary conditions in all three coordinate directions, 2)

Thin film configuration with a nonperodic boudnary in the film thickness direction

and periodic boundary conditions in the two transverse coordinate directions. For the




                                               73
bulk configuration, the simulation volume is cubic, with an edge length of approxi-

mately 11 nm at 300 K. The thin film configuration is created with a fixed, square

base measuring 8 nm by 8 nm and a film thickness that measures approximately 21

nm at 300 K.

5.2.3    Creation of Binary Silica Substrate - Supported Polystyrene
         Thin Film System

   To create the simulation system representing the supported polystyrene thin film,

the film created as described in the previous section is interfaced with the silica

slab created as described in Section 5.2.1.     To this end, the simulation volume

containing the equilibrated film is translated in the thickness direction a distance

d=hf ilm /2+hslab /2+0.4 nm and combined with the simulation volume containing the

silica slab. The combined system is then re-equilibrated for an additional 100 ps at

500 K to create a relaxed polymer-substrate interface.

5.2.4    Simulation Protocol for Determination of Polystyrene
         Tg

   The bulk PS, freestanding PS film, and supported PS film configurations are all

subjected to similar simulation protocols for cooling to measure Tg . The polymer

samples, equilibrated at 500 K, are continuously cooled to 250 K over simulation

times ranging between 1 and 24 ns. This range of simulation times corresponds to

cooling rates between 0.01 and 0.2 K/ps. From these simulations, the effects of cooling

rate, geometric configuration, and substrate interaction are analyzed individually,

as described in the following section. For the bulk PS system, periodic boundary

conditions are employed in all three coordinate directions, and the isobaric-isothermal



                                          74
ensemble is used. For the thin film configurations, a nonperiodic boundary condition

is employed in the direction of the film thickness, while periodic boundary conditions

are used for the two transverse directions. The NVT ensemble is used.


5.3    Simulation Results: Mobility-Based Polystyrene Tg Mea-
       surements

   To determine Tg , position data is collected for each united atom 30 ps during

the simulation runs described in the previous section. From this output, the mean

squared displacement is calculated for t=300 ps via Equation 3.4. This type of data

has been used successfully in such recent works as those of Doi et. al.[105] and

Baljon et. al.[106] to measure Tg as the intersection of two lines fit to the high

and low temperature MSD data. For the longest 24 ns simulations, this results in

the measurements presented in Figure 5.3. Since the value of 414 K for bulk Tg is


                                             0.4
                                                           Free      Supported   Bulk
                     MSD(t = 300ps) [nm2 ]




                                             0.3
                                                    Tg(free) =386


                                             0.2    T (supp) =396
                                                     g


                                                    Tg(bulk) =414
                                             0.1



                                              0
                                              300        350           400              450
                                                               Temperature [K]




Figure 5.3: Comparison of Tg calculations for bulk, freestanding thin film and sup-
ported thin film configurations.




                                                                    75
well above the experimenally measured 371 K, the results of this simulation must be

compared with the results of other simulations at different cooling rates in order to

extrapolate the Tg values to experimentally accessible cooling rates.

   As suggested in the work of Buchholz and co-workers,[107] a Vogel-Fulcher-Tammann[108,

109] functional form is assumed for the dependence of the relaxation time τ on the

temperature T :
                                                   EA
                              τ (T ) = τ∞ exp                                     (5.1)
                                                (T − T0 )

Here, τ∞ and EA are constants and T0 is the Vogel-Fulcher-Tammann temperature.

This, in turn, can be used to derive an expression for the glass trasition temperature

in terms of the cooling rate γ as

                                                   B
                                Tg (γ) = T0 −                                     (5.2)
                                                ln (Aγ)

where T0 is a fit parameter corresponding to the glass transition temperature in

the limit of slow cooling, and A and B are empirically determined.[110] For the

freestanding and supported thin films in this study, the values of the parameters

are determined using a nonlinear least sqaures fit to Equation 5.2, yielding the result

shown in Figure 5.4. The fitted values obtained for T0 is taken as an estimate of Tg for

PS samples cooled at an experimentally accessible rate. The value of 356K estimated

for Tg of a supported film of similar dimensions and material is in agreement with

the experimental measurments in the literature.[111, 112, 113]Some confidence in the

accuracy of this estimate is suggested by the estimate of 375 K given by the model

for bulk polystyrene, which is close to the range of 370-375 K reported experimental

values.[111, 110, 112]. However, unlike the assertion made in the paper by Keddie,

Jones, and Cory, it is shown here that the interactions between the polystyrene thin

                                          76
                                           Freestanding Film
                                           Supported Film
                                    460    Bulk


                                    440




                      T g (γ) (K)
                                    420


                                    400


                                    380
                                               −2                       −1
                                          10                           10
                                                           γ (K/p s)




Figure 5.4: Measured values of Tg for the freestanding and supported films across a
number of cooling rates. The dashed lines are fit to the data using Equation 5.2. The
values of the coefficients are A=1.05±0.06 ps/K and B=181.0±12.6 K. The values of
T0 are 342 K for the freestanding film, 356 K for the supported film, and 375 for the
bulk configuration.




film and the silica substrate have a significant impact on Tg . This results in a 14 K

difference in the estimated Tg for the freestanding and supported thin films.

   A possible reason for this observed effect is that the atomic mobility in the sup-

ported film is significantly reduced when compared to similar measurements made

in the freestanding film. This is shown in Figure 5.5. Analysis of the data in Fig-

ure 5.5(b) clearly shows that the MSD is lower in the polymer surface region near the

silica substrate than in the polymer free surface region. In fact, in the region within

1 nm of the silica substrate, the mobility of the atoms is reduced by more than an

order of magnitude. The polymer atoms in this region essentially behave as a glass

at temperatures upt to 500 K, well above Tg . Previously, it was assumed[111] that

this reduction in mobility was not significant enough to affect the measurement of Tg

for the entire film. This is a reasonable assumption, epecially given the assumption

                                                           77
                           0.75                                                                       0.75
                                                  Surface   Core                                                      Surface    Core        Interface




   MSD(t = 300ps) [nm2 ]




                                                                              MSD(t = 300ps) [nm2 ]
                                                                                                                Ts =387
                                                                                                                 g
                            0.5                                                                        0.5
                                         s
                                        T =366
                                         g
                                                                                                                Tc =392
                                                                                                                 g

                                        Tc =390
                                         g
                           0.25                                                                       0.25      Ti =424
                                                                                                                 g




                             0                                                                          0
                                  300        350        400        450                                  300          350           400               450
                                              Temperature [K]                                                              Temperature [K]

                                   (a) Freestanding thin film                                                  (b) Supported thin film



Figure 5.5: Comparison of mobility-based Tg calculations for the region near the
polymer free surface and for the core region in a freestanding polystyrene film, and
similar analyses for the free surface, interface, and core regions for a supported film
of identical composition.




of a limited extent for the region of surface-impacted mobility. However, this study

shows that the interaction of the polymer and substrate near the interface signifi-

cantly impacts Tg , as shown in Figure 5.5. Thus, a careful analysis of the mobility

characteristics is required to interpret these results.

     In Figure 5.6(a), it is shown that the extent of the free surface-affected region

of enhanced mobility approaches 5 nm in thickness at elevated temperatures, and

persists even at temperatures below Tg , as seen in Figure 5.6(b). Similarly, in the

supported film configuration, a region of reduced mobility extends from the polymer-

substrate interface approximately 4-6 nm into the film.




                                                                         78
                           0.4                                                                             0.2
                                             Freestanding   Supported                                                        Freestanding   Supported




   MSD(t = 300ps) [nm2 ]
                           0.3




                                                                                  MSD(t = 300ps) [nm2 ]
                                                                                                          0.15



                           0.2                                                                             0.1



                           0.1                                                                            0.05



                            0                                                                               0
                                 −10    −5           0           5      10                                       −10    −5            0          5      10
                                             Layer Position [nm]                                                             Layer Position [nm]

                                       (a) MSD at 400 K                                                                (b) MSD at 350 K



Figure 5.6: Layerwise mobility for freestanding and supported polystyrene films at
two temperatures. The solid horizontal line corresponds to the MSD of the bulk
system at the given temperature. For the supported film, the position of ∼ -11 nm
corresponds to the surface of the silica slab substrate.




5.4                              Conclusions

     For films with a thickness of 20 nm, the free surface and interface regions in which

the mobility is obviously affected comprises up to 40% of the total material in the

film. This distance is 3-5 times the distance over which interactions with the silica

substrate are calculated, indicating that an explanation for the reduced mobility must

go beyond simply considering the local energetic contributions from the polymer-

substrate interactions. Additionally, the relatively constant value of the MSD in the

core of the film is higher than the bulk, as seen in Figure 5.6. This suggests that the

existence of a free surface impacts the mobility of atoms throughout the entire film.

Unlike the observed increase in mobility in the film core, the relatively stable density

value observed in the film core is not significantly different than the bulk density. As

shown in Figure 4.4 in Chapter 4 and in Figure 5.7, the polymer density in the core

                                                                             79
                                             1.5




                      Layer Density (g/cc)
                                              1




                                             0.5




                                              0
                                                   −10   −5           0           5   10
                                                              Layer Position (nm)




Figure 5.7: Layerwise polymer densities for the freestanding and supported thin film
systems. The horizontal solid line represents the density measured in the bulk system.




of the film is comparable to the value observed in the bulk system.

   For the freestanding thin film, the free-surface enhanced mobility contributes to

a 33 K depression in Tg from the 375 K value for the bulk configuration to 342 K for

the film. In the case of the supported thin film, competition between the regions of

enhanced mobility at the free surface and reduced mobility at the interface surface

results in a measured Tg value of 356 K. This is 19 K lower than the bulk value, but

14 K higher than the Tg of the freestanding film.

   Thus, the Tg of the systems considered–freestanding thin film, supported thin

film, and bulk configurations–are ordered as Tf ree < Tsupported < Tbulk . This ordering
                                            g        g            g

corresponds to the inverse of the relative ordering of the average atomic mobilities

of these systems, where MSDf ree > MSDsupported > MSDbulk . This is consistent

with the hypothesis that Tg is the point at which the characteristic time is equal to

the structural relaxation time τ , given by Equation 5.1, suggesting that the overall



                                                                  80
enhanced mobility of the atoms in a thin film configuration is responsible for the

observed depression in Tg . The fact that a supported thin film exhibits a greater

reduction in Tg than the observed reduction for an identical freestanding film–as noted

by Forrest and Dalnoki-Verres[114]–is thus also explained in terms of the competing

dynamics near the substrate and free surface.

   There is an apparent contradiction between the results of this study and some

previously published results. All the results agree in that there is a reduction in the

overall film Tg . However, some observations suggest that the polymer near the sub-

strate undergoes a reduction in Tg similar to that observed near the free surface.[115]

This seems to contradict the finding in this study, that the Tg for a 1 nm thick layer

near the substrate is actually higher than that observed in the bulk material. For

example, Ellison and Torkelson[115] have found that the Tg measured in a 12-14 nm

thick layer near the substrate exhibits a 0-11 K reduction in Tg from the bulk value.

This result is obtained with 24-280 nm thick polystyrene films where the molecular

weight of the polymer is in the range of 200-500 kg/mol. For the thinnest film (24

nm) in that study, an 11 K reduction is seen in two 12 nm layers. For a slightly

thicker film (36 nm), the 14 K reduction in Tg for the 12 nm layer at the free surface

is greater than the 4 K reduction observed in the 12 nm layer near the substrate.

Similar values for the Tg for the film in this study calculated via Equation 5.2 show

that the polymer exhibits a 24 K reduction in Tg for a 10 nm thick layer near the

free surface, and a smaller 15 K reduction in Tg for a 10 nm thick layer near the

substrate. These results agree to an appropriate degree with those of Ellison and

Torkelson, given the differences in molecular weight and film height. The resolution

available for the analysis of the MD results allows Tg to be measured locally within 1


                                          81
nm thick layers. As seen in Figure 5.8, the distribution of Tg s is qualitatively similar


                                390


                                380


                                370




                      T g [K]
                                360


                                350


                                340
                                   0   5          10         15       20
                                       Distance From Substrate [nm]




        Figure 5.8: Layerwise Tg estimates for a 21 nm supported thin film.




to the distribution of mobilities pictured in Figure 5.6, as expected.




                                                82
          Appendix A: SYMPLECTIC INTEGRATORS




   To simulate the time-evolution of systems governed by the equations given in the

previous section, a symplectic algorithm[116] is required to conserve the Hamiltonian

of the system. For the purposes of the following development, it is assumed that all

masses are constant, and that the functional form of the potential does not vary with

time.

   Hamilton’s equations can be written as

                                  ˙
                                  q    0 −I
                                    =−      ∇H                                  (A.1)
                                  ˙
                                  p    I 0
               ∂
               ∂q                                            0 −I
where ∇ =      ∂    and (q, p) ∈ R6N . The matrix J =                is known as the
               ∂p
                                                             I 0
symplectic identity.[26] Letting x = (p, q), a compact form of Hamilton’s equations is


                                    x = J −1 ∇H(x, t)
                                    ˙                                           (A.2)


A symplectic map, then, is a transformation φ : R6N → R6N , represented by a

6N × 6N matrix F , with x = φ(X) = F X, where the matrix F satisfies


                                     F T J −1 F = J −1                          (A.3)


so that det (F ) = 1, i.e., Liouville’s theorem holds.[26]



                                            83
   The Velocity Verlet (VV) algorithm[29] is a second-order integrator given by

        ∆t             ∆t
 ˜
 p t+           =P −      ∇U (Q)
        2              2
                                      ∆t                      (∆t)2 −1
               q = Q + ∆tM −1 p t +        = Q + ∆tM −1 P −        M ∇U (Q) (A.4)
                                      2                         2
                         ∆t       ∆t              (∆t)2          ∆t
             p=p t+           −      ∇U (q) = P −       ∇U (Q) −    ∇U (q)
                         2        2                 2            2

                                                  ˜
where Q and P are the momentum vectors at time t, p is the half-step momenta, q

and p are the position and momentum vectors at time t + ∆t and M −1 is the inverse

of the 3N × 3N mass matrix.

   The partial derivatives that form the matrix F are given by

     ∂q           (∆t)2 −1
          =I−          M ∇ ⊗ U (Q)
     ∂Q             2
     ∂q
          = ∆tM −1
     ∂P
                                                                                    (A.5)
     ∂p      ∆t                     (∆t)3 −1
          = − ∇ ⊗ [U (Q) + U (q)] +      M [∇ ⊗ U (Q)] [∇ ⊗ U (q)]
     ∂Q      2                        4
     ∂p        (∆t)2 −1
          =I−       M [∇ ⊗ U (Q)]
     ∂P          2
where I is the 3N × 3N identity matrix, so that
                                                                               
                         2

                I − (∆t) M −1 ∇ ⊗ U (Q)
                       2
                                                          ∆tM     −1
                                                                                
  F =           − ∆t ∇ ⊗ [U (Q) + U (q)]
                                                                                   (A.6)
                   2                                 (∆t)2
                                                                               
                 3                             I−          M −1   [∇ ⊗ U (Q)]
           + (∆t) M −1 [∇ ⊗ U (Q)] [∇ ⊗ U (q)]
               4
                                                       2


which can easily be shown to satisfy Equation A.3. Thus, the Verlet algorithm is

symplectic, and as such is an acceptable choice for the integration required in MD

simulations.




                                           84
     Appendix B: DETERMINATION OF A MAXIMUM
  TIMESTEP FOR STABILITY OF THE INTEGRATION
         ALGORITHM IN MOLECULAR DYNAMICS
                               SIMULATIONS




   The use of Molecular Dynamics (MD) to study time-dependent phenomena cur-

rently exhibits a computational cost that prohibits a full description of full-blown

devices, and limits study to atomic-scale regions of interest. Furthermore, study is

limited to timescales orders of magnitude below those on which phenomena of interest

occur. Given that there is a need to study a wide range of dynamic material response

at the device level, it is imperative to develop methods that allow vast reductions in

the computational load.

   Within the MD framework, there are two broad categories that comprise the

bulk of attempts to provide the required advancements in scale: 1) Reduction of the

number of degrees of freedom in the system; 2) Increase in the allowable timestep for

integration. One particular method that combines these two approaches is to ’coarse-

grain’ explicit atoms into ’united-’ or ’super-’ atoms. This technique reduces the

total degrees of freedom for a system by essentially removing a number of atoms and

has the potential to increase the allowable timestep by suppressing higher frequency

oscillations in the system.

                                         85
   For the explicit atom representation of polystyrene, the C-H bond is the highest

frequency oscillator, and the maximum allowable timestep is determined by the sta-

bility of the computation of the dynamics of this oscillation. For the united atom

system, the CH-CH bond becomes the highest frequency oscillator, so the allowable

timestep is limited by that frequency. To determine the acceptable range of timesteps,

it is sufficient to consider the eigenvalues of the matrix F given by Equation A.6 in

Appendix A. The eigenvalues can be determined by solving the equation


                                    det (F − λI) = 0                                (B.1)


for the eigenvalues λi . For stable integration of the dynamics of a single harmonic

bond, the eigenvalues must be a complex conjugate pair.[117] Substituting Equa-

tion A.6 into Equation B.1 gives

                               2(∆t)2
                            λ +       k−2 λ+1= 0                                    (B.2)
                                  µ

                                                       m1 m2
where k is the bond energy paramter, and µ =           m1 +m2
                                                              ,   where m1 and m2 are the

masses of the two bonded atoms. For λ1 , λ2 to form a complex conjugate pair,
                                               2
                                   (∆t)2
                                         k−2       −4<0                             (B.3)
                                     µ

must hold. This implies that
                                                   µ
                                      ∆t < 2                                        (B.4)
                                                   k

Substituting the actual values for the C-H and CH-CH bonds from Tables 2.1 and

2.2, the maximum allowable timestep for the CH-CH bond is 10.87 fs, or 2.06 times

the maximum allowable timestep of 5.28 fs for the C-H bond.




                                          86
                                Bibliography



 [1] M. P. Allen and D. J. Tildesley. Computer simulation of liquids. Oxford Uni-
     versity Press, 1987.

 [2] Y.C. Jean. Positron annihilation spectroscopy for chemical analysis: A
     novel probe for microstructural analysis of polymers. Microchemical Journal,
     42(1):72–102, 1990.

 [3] R. G. Wissinger and M. E. Paulaitis. Molecular thermodynamic model for sorp-
     tion and swelling in glassy polymer-co2 systems at elevated pressures. Industrial
     and Engineering Chemistry Research, 30:842, 1991.

 [4] Hongmin Chen, Mei-Ling Cheng, Y. C. Jean, L. James Lee, and Jintao Yang.
     Effect of co2 exposure on free volumes in polystyrene studied by positron an-
     nihilation spectroscopy. Journal of Polymer Science Part B: Polymer Physics,
     46:388–405, 2008.

               e
 [5] Suichi Nos´. A unified formulation of the constant temperature molecular dy-
     namics methods. Journal of Chemical Physics, 81(1):511–518, 1984.

 [6] I. T. Todorov and W. Smith. Dl poly 3: the ccp5 national uk code for molecular-
     dynamics simulations. Proceedings of the Royal Society of London A, 362:1835–
     1852, 2004.

 [7] Sandia National Laboratories. LAMMPS Users Manual, January 2010.

 [8] A. Srivastava and S. Ghosh. Experimentally validated md model for glass
     transition temperature and free volume fraction of polystyrene. International
     Journal for Multiscale Computational Engineering, 8(6), 2010.

 [9] Ali A. Hassanali and Sherwin J. Singer. Model for the water-amorphous silica
     interface: The undissociated surface. Journal of Physical Chemistry Part B,
     111(38), 2007.

[10] J. Han and R. H. Boyd. Molecular packing and small penetrant diffusion in
     polystyrene: A molecular dynamics simulation study. Polymer, 37(10):1797–
     1804, 1996.

                                         87
[11] George Kaminski and William L. Jorgensen. Performance of the amber94,
     mmff94, and opls-aa force fields for modeling organic liquids. Journal of Physical
     Chemistry, 100:18010–18013, 1996.

[12] Jean-Paul Ryckaert, Giovanni Ciccotti, and Herman J.C. Berendsen. Numerical
     integration of the cartesian equations of motion of a system with constraints:
     Molecular dynamics of n-alkanes. Journal of Computational Physics, 23:327–
     341, 1977.

[13] Wilfred F. van Gunsteren and Martin Karplus. Effect of constraints on the
     dynamics of macromolecules. Macromolecules, 15:1528–1544, 1982.

[14] S.Toxvaerd. Molecular dynamics calculation of the equation of state of alkanes.
     Journal of Chemical Physics, 6:4290, 1990.

[15] J. Han, R.H. Gee, and R.H. Boyd. Glass transition temperatures of polymers
     from molecular dynamics simulations. Macromolecules, 27:7781–7784, 1994.

[16] Collin D. Wick, Marcus G. Martin, and J. Ilja Siepmann. Transferable po-
     tentials for phase equilibria. 4. united-atom description of linear and branched
     alkenes and alkylbenzenes. J. Phys. Chem. B, 104:8008–8016, 2000.

[17] V. A. Harmandaris, N. P. Adhikari, N. F. A. van der Vegt, and K. Kremer. Hier-
     archical modeling of polystyrene: From atomistic to coarse-grained simulations.
     Macromolecules, 39:6708–6719, 2006.

[18] J. G. Harris and K. H. Yung. Carbon dioxide’s liquid-vapor coexistence curve
     and critical properties as predicted by a simple molecular model. Journal of
     Physical Chemistry, 99:12021–12024, 1995.

[19] Moumita Saharay and Sundram Balasubramanian. Ab initio molecular-
     dynamics study of supercritical carbon dioxide. Journal of Chemical Physics,
     120(20):9694, 2004.

[20] Patrick Etesse, James A. Zega, and Riki Kobayashi. High pressure nuclear
     magnetic resonance measurement of spin-lattice relaxation and self-diffusion in
     carbon dioxide. Journal of Chemical Physics, 97(3):2022–2029, 1992.

[21] Lester Guttman and Shafiqur M. Rahman. Simulation of the structure of amor-
     phous silicon dioxide. Physical Review B, 37(5):2657–2668, 1988.

[22] B. W. H. van Beest, G. J. Kramer, and R. A. van Santen. Force fields for silicas
     and aluminophosphates based on ab initio calculations. Physical Review Lettes,
     64(16):1955–1958, 1990.



                                        88
[23] P. N. Keating. Effect of invariance requirements on the elastic strain en-
     ergy of crystals with application to the diamond structure. Physical Review,
     145(2):637–645, 1966.

[24] B. P. Feuston and S. H. Garofalini. Empirical three-body potential for vitreous
     silica. Journal Of Chemical Physics, 89(9):5818–5824, 1988.

[25] Z. Jiang and R. A. Brown. Modelling oxygen defects in silicon crystals using an
     empirical interatomic potential. Chemical Engineering Science, 49(17):2991–
     3000, 1994.

[26] Sherwin J. Singer. Lecture notes, chemistry 996, the ohio state university. 2010.

[27] Maurice L. Huggins and Joseph E. Mayer. Interatomic distances in crystals of
     the alkali halides. Journal Of Chemical Physics, 1:643–646, 1933.

[28] K. Vollmayr-Lee, J. A. Roman, and J. Horbach. Aging to equilibrium dynamics
     of sio2 . Physical Review E, 81:061203, 2010.

[29] I. T. Todorov and W. Smith. The DL POLY 3 User Manual. Daresbury Lab-
     oratory, September 2007.

[30] J.E. Lennard-Jones. On the determination of molecular fields.-ii. from the equa-
     tion of state of a gas. Proceedings of the Royal Society of London, Series A,
     106:463, 1924.

[31] T. A. Halgren. Representation of van der waals (vdw) interactions in molecular
     mechanics force fields: Potential form, combination rules, and vdw parameters.
     Journal of the American Chemical Society, 114:7827–7843, 1992.

[32] Isaac C. Sanchez and Robert H. Lacombe. Statistical thermodynamics of poly-
     mer solutions. Macromolecules, 11(6):1145–1156, 1978.

[33] Chongli Zhong and Hirokatsu Masuoka. Modeling of gas solubilities in polymers
     with cubic equation of state. Fluid Phase Equilibra, 144:49–57, 1998.

[34] Yun Kyung Shin, Ali A. Hassanali, and Sherwin J. Singer. Biomolecules at
     the amorphous silica/water interface: model development and application to
     binding and fluorescence anistropy of peptides. 2010.

[35] P. Flory. Statistical Mechanics of Chain Molecules. Interscience, 1969.

[36] R.F. Rapold and U.W. Suter. Conformational characteristics of polystyrene.
     Macromolecular Theory and Simulations, 3:1–17, 1994.



                                         89
[37] Roland F. Rapold, Ulrich W. Suter, and Doros N. Theodorou. Static atomistic
     modelling of the structure and ring dynamics of bulk amorphous polystyrene.
     Macromolecular Theory and Simulations, 3:19–43, 1994.

[38] Jonathan I. McKechnie, David Brown, and Julian H. R. Clarke. Methods of
     generating dense relaxed amorphous polymer samples for use in dynamic sim-
     ulations. Macromolecules, 25:1562–1567, 1992.

[39] David Brown, Julian H. R. Clarke, Motoi Okuda, and Takao Yamazaki. The
     preparation of polymer melt samples for computer simulation studies. Journal
     of Chemical Physics, 8:6011–6018, 1994.

[40] R. Auhl, R. Everaers, G.S. Grest, K. Kremer, and S.J. Plimpton. Equilibration
     of long chain polymer melts in computer simulations. Journal of Chemical
     Physics, 119(24):12718–12728, 2003.

[41] Andrew I. Cooper. Polymer synthesis and processing using supercritical carbon
     dioxide. Journal of Materials Chemistry, 10:207–234, 2000.

[42] S. Nalwade, F. Picchioni, and L. Janssen. Supercritical carbon dioxide as a
     green solvent for processing polymer melts: Processing aspects and applications.
     Progress in Polymer Science, 31:19, 2006.

[43] Eizo Sada, Hidehiro Kumazawa, Hiroshi Yakushiji, Yukio Bamba, Kazuhiko
     Sakata, and Shao-Ting Wang. Sorption and diffusion of gases in glassy polymers.
     Industrial and Engineering Chemistry Research, 26:433–438, 1987.

[44] S.G.Kazarian, M.F.Vincent, C.L.Liotta, and C.A.Eckert. Specific intermolec-
     ular interaction of carbon dioxide with polymers. J. Am. Chem. Soc., 9:1729,
     1996.

[45] Y. Sato, M. Yurugi, K. Fujiwara, S. Takishima, and H. Masuoka. Solubili-
     ties of carbon dioxide and nitrogen in polystyrene under high temperature and
     pressure. Fluid Phase Equilibra, 125:129–138, 1996.

[46] C. J. Peng, H. L. Liu, and Y. Hu. Gas solubilities in molten polymers based on
     an equation of state. Chemical Engineering Science, 56:6967–6975, 2001.

[47] Z. Dong and J. R. Fried. Statistical thermodynamics of the glass transition:
     1. effect of pressure and diluent concentration. Computational and Theoretical
     Polymer Science, 7(1):53–64, 1997.

[48] Y. Yang. Carbon dioxide assisted polymer micro/nanofabrication. PhD thesis,
     The Ohio State University, 2005.



                                        90
[49] Y.Yang, D.Lui, Y.Xie, L.J.Lee, and D.L.Tomasko. Low temperature fusion of
     polymeric nanostructures using carbon dioxide. Advanced Materials, 19:251–
     254, 2007.

[50] H. Abramowitz, P.S. Shah, P.F. Green, and K.P. Johnston. Welding colloidal
     crystals with carbon dioxide. Macromolecules, 37:7316, 2004.

[51] R. G. Wissinger and M. E. Paulaitis. Glass transitions in polymer / co2 mixtures
     at elevated pressures. Journal of Polymer Science Part B: Polymer Physics,
     29:631–633, 1991.

[52] P.D. Condo, I.C. Sanchez, C.G. Panayiotou, and K.P. Johnston. Glass tran-
     sition behavior including retrograde vitrification of polymers with compressed
     fluid diluents. Macromolecules, 25:6119, 1992.

[53] P.D.Condo, D.R.Paul, and K.P.Johnston. Glass transitions of polymers with
     compressed fluid diluents: Type i1 and i11 behavior. Macromolecules, 27:365,
     1994.

[54] C. Kwag, C.W. Manke, and E. Gulari. Effects of dissolved gas on viscoelastic
     scaling and glass transition temperature of polystyrene melts. Industrial and
     Engineering Chemistry Research, 40:3048, 2001.

[55] Hisao Takeuchi and Keiji Okazaki. Molecular dynamics simulation of diffusion
     of simple gas molecules in a short chain polymer. Journal of Chemical Physics,
     92(9):5643–5652, 1990.

[56] Dumitru Pavel and Robert Shanks. Molecular dynamics simulation of diffusion
     of o2 and co2 in amorphous poly(ethylene terephthalate) and related aromatic
     polyesters. Polymer, 44(21):6713–6724, 2003.

                                    u
[57] Hossein Eslami and Florian M¨ ller-Plathe. Molecular dynamics simulation of
     sorption of gases in polystyrene. Macromolecules, 40:6413–6421, 2007.

[58] Y.C. Jean, T.C. Sandreczki, and D.P. Ames. Positronium annihilation in amine-
     cured epoxy polymers. Journal of Polymer Science Part B, 24:1247, 1986.

[59] H. Nakanishi, S. Wang, and Y. C. Jean. Positron Annihilation Studies of Fluids,
     page 292. World Science: Singapore, 1988.

[60] X. Hong, Y.C. Jean, Hsinjin Yang, S.S. Jordon, and W.J. Koros. Free-volume
     hole properties of gas-exposed polycarbonate studied by positron annihilation
     lifetime spectroscopy. Macromolecules, 29:7859, 1996.




                                        91
[61] J. Yuan, H. Cao, E.W. Hellmuth, and Y.C. Jean. Subnanometer hole proper-
     ties of co2 -exposed polysulfone studied by positron annihilation lifetime spec-
     troscopy. Journal of Polymer Science Part B, 36:3049, 1998.

[62] S.E. Baltazar, A.H. Romero, J.L. Rodriguez-Lopez, H. Terrones, and R. Mar-
     tonak. Assessment of isobaric-isothermal (npt) simulations for finite systems.
     Computational Material Science, 37:526, 2006.

[63] Richard Vawter. Physicsnet - learning physics on the web. http://faculty.
     wwu.edu/vawter/PhysicsNet/Topics/Thermal/vdWaalEquatOfState.html,
     September 2010.

[64] Jeffrey Ellis, David L. Tomasko, and Fariba Dehghani. Novel dense co2
     technique for β-galactosidase immobilization in polystyrene microchannels.
     Biomacromolecules, 9:1027–1034, 2008.

[65] D. Gersappe and M. O. Robbins. Where do polymer adhesives fail? Europhysics
     Letters, 48(2):150–155, 1999.

[66] R.G. Wissinger and M.E. Paulaitis. Swelling and sorption in polymer-co2
     mixtures at elevated pressures. Journal of Polymer Science Part B: Polymer
     Physics, 25(12):2497–2510, 1987.

[67] J. S. Chiou, J. W. Barlow, and D. R. Paul. Plasticization of glassy polymers
     by co2 . Journal of Applied Polymer Science, 30(6):2633–2642, 1985.

[68] A. Srivastava, C. Alleman, S. Ghosh, and L. J. Lee. Molecular dynamics sim-
     ulation based evaluation of glass transition temperatures of polystyrene in the
     presence of carbon dioxide. Modelling and Simulation in Materials Science and
     Engineering, 18(6):065003, 2010.

[69] Theodora Spyriouni, Georgios Boulougouris, and Doros N. Theodorou. Pre-
     diction of sorption of co2 in glassy atactic polystyrene at elevated pressures
     through a new computational scheme. Macromolecules, 42:1759–1769, 2009.

[70] Y.C. Jean, Renwu Zhang, H. Cao, Jen-Pwu Yuan, Chia-Ming Huang,
     B. Nielsen, and P. Asoka-Kumar. Glass transition of polystyrene near the sur-
     face studied by slow-positron-annihilation spectroscopy. Physical Review B,
     56(14):R8459–R8462, 1997.

[71] Noriaki Satomi, Keiji Tanaka, Atsushi Takahara, and Tisato Kajiyama. Surface
     molecular motion of monodisperse α,ω-diamino-terminated and α,ω-dicarboxy-
     terminated polystyrenes. Macromolecules, 34:8761–8767, 2001.



                                        92
[72] Rajaram A. Pai, Raashina Humayun, Michelle T. Schulberg, Archita Sengupta,
     Jia-Ning Sun, and James J. Watkins. Mesoporous silicates prepared using pre-
     organized templates in supercritical fluids. Science, 303:507, 2004.

[73] Young Hwa Kim and Richard P. Wool. A theory of healing at a polymer-polymer
     interface. Macromolecules, 16:1115–1120, 1983.

[74] Ralf Schnell, Manfred Stamm, and Constantino Creton. Direct correlation
     between interfacial width and adhesion in glassy polymers. Macromolecules,
     31:2284–2292, 1998.

[75] Kei-Ichi Akabori, Daisuke Baba, Kazuhiro Koguchi, Keiji Tanaka, and Toshi-
     hiko Nagamura. Relation between the adhesion strength and interfacial width
     for symmetric polystyrene bilayers. Journal of Polymer Science Part B: Poly-
     mer Physics, 44:3598–3604, 2006.

[76] Manish Talreja, Isamu Kusaka, and David L. Tomasko. Density functional
     approach for modeling co2 pressurized polymer thin films in equilibrium. The
     Journal of Chemical Physics, 130(8):084902, 2009.

[77] Yuri M. Boiko and Robert E. Prud’homme. Strength development at the inter-
     face of amorphous polymers and their miscible blends, below the glass transition
     temperature. Macromolecules, 31:6620–6626, 1998.

[78] Yuri M. Boiko, Anders Bach, and Jorgen Lyngaae-Jorgensen. Self-bonding in
     an amorphous polymer below the glass transition: A t-peel test investigation.
     Polymer Science Part B: Polymer Physics, 42:1861–1867, 2004.

[79] Xiaomin Zhang, Shigeru Tasaka, and Norihiro Inagaki. Surface mechanical
     properties of low-molecular-weight polystyrene below its glass-transition tem-
     peratures. Journal of Polymer Science Part B: Polymer Physics, 38:654–658,
     2000.

[80] Jintao Yang, Chuntai Liu, Yong Yang, Bin Zhu, L. James Lee, Hongmin Chen,
     and Y. C. Jean. Analysis of polystyrene surface properties on thin film bonding
     under carbon dioxide pressure using nanoparticle embedding technique. Poly-
     mer Science Part B: Polymer Physics, 47:1535–1542, 2009.

[81] R. P. Wool. Polymer Interfaces: Structure and Strength. Hanser, 1995.

[82] Daisuke Kawaguchi, Keiji Tanaka, and Tisato Kajiyama. Mobility gradient
     in surface region of monodisperse polystyrene films. Macromolecules, 36:1235–
     1240, 2003.

[83] Shin Kawana and Richard A. L. Jones. Character of the glass transition in thin
     supported polymer films. Physical Review E, 63:021501, 2001.

                                        93
[84] Stephen Prager and Matthew Tirrell. The healing process at polymer-polymer
     interfaces. Journal of Chemical Physics, 75(10):5194–5198, 1981.

[85] Jr. Prince E. Rouse. A theory of the linear viscoelastic properties of dilute
     solutions of coiling polymers. The Journal of Chemical Physics, 21(7):1272–
     1280, 1953.

[86] D. H. Tsai. The virial theorem and stress calculation in molecular dynamics.
     Journal of Chemical Physics, 70(03):1375–1382, 1979.

[87] Nikhil Chandra Admal and E. B. Tadmor. A unified interpretation of stress in
     molecular systems. Journal of Elasticity, 100:63–143, 2010.

[88] E. Raphael and P. G. de Gennes. Rubber-rubber adhesion with connector
     molecules. Journal of Physical Chemistry, 96(10):4002–4007, 1992.

[89] Hong Ji and P. G. de Gennes. Adhesion via connector molecules: The many-
     stitch problem. Macromolecules, 26:520–525, 1993.

             e
[90] A. Jim´nez, G. L. Thompson, M. A. Matthews, T.A. Davis, K. Crocker, J.S.
     Lyons, and A. Trapotsis. Compatibility of medical grade polymers with dense
     co2 . Journal of Supercritical Fluids, 42:366–372, 2007.

[91] Wen-Chou V. Wang, Edward J. Kramer, and Wolfgang H. Sachse. Effects of
     high-pressure co2 on the glass transition temperature and mechanical prop-
     erties of polystyrene. Journal of Polymer Science: Polymer Physics Edition,
     20(8):1371–1384, 1982.

[92] Jessica M. Torres, Christopher M. Stafford, and Bryan D. Vogt. Elastic modulus
     of amorphous polymer thin films: Relationship to the glass transition temper-
     ature. ACS Nano, 3(9):2677–2685, 2009.

[93] Francis W. Starr, Thomas B. Schrøder, and Sharon C. Glotzer. Interface and
     surface effects on the glass-transition temperature in thin polymer films. Fara-
     day Discussions, 98:219–230, 1994.

[94] Mehrzad Mortezaei, Gholamali Farzi, Mohammad Reza Kalaee, and Mahmood
     Zabihpoor. Evaluation of interfacial layer properties in the polystyrene/silica
     nanocomposite. Journal of Applied Polymer Science, 119:2039–2047, 2011.

[95] Perla Rittigstein and John M. Torkelson. Polymernanoparticle interfacial in-
     teractions in polymer nanocomposites: Confinement effects on glass transition
     temperature and suppression of physical aging. Journal of Polymer Science
     Part B, 44(20):29352943, 2006.



                                        94
 [96] Joseph L. Keddie, Richard A. L. Jones, and Rachel A. Cory. Interface and sur-
      face effects on the glass-transition temperature in thin polymer films. Faraday
      Discussions, 98:219–230, 1994.

 [97] J.A. Forrest. A decade of dynamics in thin films of polystyrene: Where are we
      now? The European Physical Journal E, 8:261–266, 2002.

 [98] L. V. Woodcock, C. A. Angell, and P. Cheeseman. Molecular dynamics studies
      of the vitreous state: Simple ionic systems and silica. Journal Of Chemical
      Physics, 65(4):1565–1577, 1976.

 [99] T. F. Soules. A molecular dynamic calculation of the structure of sodium silicate
      glasses. Journal Of Chemical Physics, 71(11):4570–4578, 1979.

[100] P. Masini and M. Bernasconi. Ab initio simulations of hydroxylation and de-
      hydroxylation reactions at surfaces: amorphous silica and brucite. Journal of
      Physics: Condensed Matter, 14:4133–4144, 2002.

[101] Mark Wilson and Tiffany R. Walsh. Hydrolysis of the amorphous silica surface.
      i. structure and dynamics of the dry surface. Journal Of Chemical Physics,
      113(20):9180–9190, 2000.

[102] Tiffany R. Walsh, Mark Wilson, and Adrian P. Sutton. Hydrolysis of the amor-
      phous silica surface. ii. calculation of activation barriers and mechanisms. Jour-
      nal Of Chemical Physics, 113(20):9191–9201, 2000.

[103] L. T. Zhuravlev. Characterization of amorphous silica surface. Reaction Kinetics
      and Catalysis Letters, 50(1-2):15–25, 1993.

[104] Ralph K. Iler. The Chemistry of Silica: Solubility, Polymerization, Colloid and
      Surface Properties, and Biochemistry. John Wiley & Sons, 1979.

[105] Hiroshi Morita, Keiji Tanaka, Tisato Kajiyama, Toshio Nishi, and Masao Doi.
      Study of the glass transition temperature of polymer surface by coarse-grained
      molecular dynamics simulation. Macromolecules, 39(18):62336237, 2006.

[106] R. C. Baljon, S. Williams, N. K. Balabaev F., Paans, D. Hudzinskyy, and
      A. V. Lyulin. Simulated glass transition in free-standing thin polystyrene film.
      Journal of Polymer Science Part B: Polymer Physics, 48:11601167, 2010.

[107] Joachim Buchholz, Wolfgang Paul, Fathollah Varnik, and Kurt Binder. Cooling
      rate dependence of the glass transition temperature of polymer melts: Molecular
      dynamics study. Journal of Chemical Physics, 117(15):7364–7372, 2002.

                                   a                            a       u
[108] Hans Vogel. Das tempetatur¨bhangigkeitsgesetz der viskosit¨t von fl¨ ssigkeiten.
      Physikalische Zeitschrift, 22:645–646, 1921.

                                          95
[109] Gordon S. Fulcher. Analysis of recent measurements of the viscosity of glasses.
      Journal of the American Ceramic Society, 8(5):339–355, 1925.

[110] Alexey V. Lyulin, Nikolaj K. Balabev, and M. A. J. Michels. Molecular-weight
      and cooling-rate dependence of simulated tg for amorphous polystyrene. Macro-
      molecules, 36:8574–8575, 2003.

[111] J. L. Keddie, R. A. Jones, and R. A. Cory. Size-dependent depression of the
      glass transition temperature in polymer films. Europhysics Letters, 27:59–64,
      1994.

[112] Christopher J. Ellison, Manish K. Mundra, and John M. Torkelson. Impacts of
      polystyrene molecular weight and modification to the repeat unit structure on
      the glass transition-nanoconfinement effect and the cooperativity length scale.
      Macromolecules, 38:1767–1778, 2005.

[113] Manish K. Mundra, Christopher J. Ellison, Ross E. Behling, and John M.
      Torkelson. Confinement, composition, and spin-coating effects on the glass tran-
      sition and stress relaxation of thin films of polystyrene and styrene-containing
      random copolymers: Sensing by intrinsic fluorescence. Polymer, 47(22):7747–
      7759, 2006.

[114] James A. Forrest and Kari Dalnoki-Veress. The glass transition in thin polymer
      films. Advances in Colloid and Interface Science, 2:167–196, 2001.

[115] Cristopher J. Ellison and John M. Torkelson. The distribution of glass-
      transition temperatures in nanoscopically confined glass formers. Nature Mate-
      rials, 94:695–700, 2003.

[116] Ronald D. Ruth. A canonical integration technique. IEEE Transactions on
      Nuclear Science, NS-30(4):2669–2671, 1983.

[117] Jong-In Choe and Byungchul Kim. Determination of proper time step for molec-
      ular dynamics simulation. Bulletin of the Korean Chemical Society, 21(4):419–
      424, 2000.




                                         96

								
To top