NANO LETTERS
Computer Simulation of the Mechanical Properties of Amorphous Polymer Nanostructures
2003 Vol. 3, No. 10 1405-1410
Kevin Van Workum* and Juan J. de Pablo UniVersity of Wisconsin-Madison, Madison, Wisconsin Received June 30, 2003
ABSTRACT The mechanical properties of polymeric nanostructures have been investigated by means of Monte Carlo and molecular dynamics simulations. A continuum mechanics model was applied to determine an apparent modulus for such structures from results of virtual deformation simulations. A recently proposed method based on strain fluctuations was also used to calculate the elasticity tensor for the structures. The effect of size on the apparent modulus of ultra-small structures was explored. Our results indicate that, for a given system at a specified temperature, the modulus of a small structure can be significantly smaller than that of the bulk material. Furthermore, in small systems the elastic constants are shown to become anisotropic.
1. Introduction. Over the past few years the theoretical and numerical study of supercooled liquids and glasses has received considerable attention. Most of the published work in this field has been concerned with bulk materials, where a number of important questions remain unsolved. Experiments on confined systems1,2 and thin polymer films,3 however, have shown that confinement can considerably alter the transition from a liquid to a glassy state. The origin of the observed departures from bulk behavior remains elusive. In addition to its inherent fundamental importance, the study of glass transition and related phenomena in confined systems has important practical applications. The semiconductor industry envisages large-scale production of sub-100 nm structures soon, which are necessary to improve the performance of current integrated circuits.4 As the production of microdevices by lithographic methods enters the nanometer length scale, it is becoming increasingly clear that many of the materials involved in their fabrication behave differently than in the bulk.5,6 Some of these materials are amorphous polymeric glasses. Our current understanding of structure-property relations in polymer glasses at the nanoscale level is severely limited, and this lack of knowledge poses a major obstacle toward rational development of novel materials and processes for fabrication of ultra-small devices. Glassy materials are often described in terms of their stiffness or rigidity. The glass transition temperature is often identified as the temperature at which the viscosity of a glassforming liquid attains a value of 1013 centipoise (cp). From a practical point of view, the mechanical properties of a glass often dictate whether a specific need or application can be met. Ironically, theoretical or numerical studies of the mechanical properties of glasses have been scarce. This work fills a gap in the literature by presenting a detailed numerical 10.1021/nl034458l CCC: $25.00 Published on Web 09/04/2003
© 2003 American Chemical Society
study of the mechanical properties of a model polymeric glass and their departure from bulk behavior in nanometer-scale structures, analogous to those encountered in lithographic nanofabrication processes. More specifically, given that the glass transition temperature of thin polymer films has been found to deviate from its bulk value,3,7,8 the fundamental question addressed in this work is whether the elastic moduli of a nanometer lengthscale polymeric structure are different from those of the bulk material. A number of recent experiments have been carried out in attempts to quantify the elastic constants of small systems.9-11 These experiments, however, face considerable challenges as the dimensions of a system become small. At the 50 nm scale and below, such experiments would benefit considerably from the insights provided by molecular simulations. Recent simulations of nanostructures8 have shown that their bending modulus can be considerably smaller than that of the bulk polymer. Those calculations, however, were conducted with a model in which interactions between monomers were of the square-well type. Furthermore, elastic constants were extracted from deformation simulations under the assumption that bending-plate equations from continuum mechanics are valid at the nanoscale. This work presents results of simulations of model structures similar to those considered by Boehme and de Pablo.8 It provides a departure from that study in that a continuous potential energy function is employed, thereby permitting use of strain-fluctuation formulas and avoiding the need for continuum-mechanics assumptions. The results provided here indicate that the elastic properties of polymeric nanoscopic structures decrease significantly and become anisotropic. The outline for this paper is as follows. In sections 2 and 3 we describe the molecular simulation techniques used to
calculate the elasticity tensor from strain fluctuations in the bulk and in the nanoscopic structures. In section 4 the details of the polymer model used and the results for the bulk system are presented. Simulation details and results for the nanoscopic structures are presented in section 5, followed by a brief summary of our findings. 2. NσT Ensemble Monte Carlo Simulations. The constant temperature and stress ensemble (NσT) is a generalization of the NPT ensemble, in which the size and shape of the simulation box are allowed to change.12 System configurations are generated in the isothermal-isostress ensemble using the Metropolis Monte Carlo method.13 One Monte Carlo move consists of a random attempt to randomly displace all the particle positions, rR,i, or a random displacement of the scaling matrix, hij. Trial configurations are accepted with the following probability: Ptrans ) min [1, exp(-∆H/kBT)]
(1)
where ∆H is the change in enthalpy associated with the trial move. The enthalpy change is given by ∆H ) ∆U + P0∆V + V0σij∆ji - NkBT log(V′/V) (2) where ∆U is the change in energy of the proposed move, P0 is the external pressure, ∆V ) V′ - V is the change in volume, V0 is the volume of the reference state, σij is the imposed stress tensor, and ij is the strain tensor. Subscripts i and j refer to the Cartesian components of the coordinate system and repeated indices imply a summation. The number of particles is denoted by N, kB is the Boltzmann constant, and T is the temperature. The strain tensor is given by 1 -1 ij ) {〈h˜ 〉-1 ik Gkl〈h〉lj - δij} 2
Cijkl )
kB T [〈 〉 - 〈ij〉〈kl〉]-1 V0 ij kl
(6)
Small fluctuations of the strain tensor are characteristic of stiff materials with large elastic constants. Large fluctuations of the strain tensor correspond to soft materials with small elastic constants. 3. Elastic Bath Method. The technique outlined above was originally formulated for bulk systems under periodic boundary conditions. An extension is proposed in this section to permit its application to small samples of arbitrary shape and to improve the convergence and accuracy of the strain fluctuation formula (eq 6). For small strains, the strain energy per unit volume of a deformed system can be written in general as 1 Eˆ s ) Cijkllkji 2
(7)
For an isotropic material the strain energy per unit volume16 is given by 1 Eˆ s ) λ2ii + µ2ij 2
(8)
(3)
If a, b, and c are the principal axes of the simulation box, then hij ) {a,b,c}ij describes its size and shape. The matrix 〈h〉ij describes the reference box and 〈h˜ 〉 -1 ij is the inverse of the transpose of 〈h〉ij. The metric tensor Gij is defined as Gij ) h˜ ikhkj. The scaling matrix, hij, relates the real particle coordinates, rR,i, to the scaled particle coordinates, sR,i, according to rR,i ) hijsR,j
where ∆hmax is the maximum allowed displacement of an element of hij and ξij is a matrix of random numbers uniformly distributed between 0 and 1. To prevent rotations of the simulation box, we impose the constraint that hij ) hji. Random displacements in rR,i are made in the usual manner. The isothermal elasticity tensor can be expressed in terms of the thermal strain fluctuations at a zero applied stress according to14,15
(4)
where λ and µ are the Lame´ coefficients. If the system of interest is immersed in an elastic bath having appropriate elastic constants, the thermal strain fluctuations of the system can be reduced (or increased) significantly by arbitrarily stiffening (or softening) the resulting, ideal composite. To do this, we include the strain energy of the elastic bath from either eq 7 or eq 8 in the Monte Carlo acceptance criteria. Equation 2 then becomes ∆H ) ∆U + P0∆V + V0σij∆ji + ∆Ebs - NkBT log(V′/V) (9)
The scaled particle coordinates describe the positions of the particles in a cubic box with sides of unit length. The real particle coordinates describe the positions of the particles in a box whose size and shape are given by hij. The matrix form of hij restricts the shape of the simulation box to a general parallelepiped; only homogeneous strains are possible in our simulations. Random displacements in hij are made in the following way: h′ij ) hij + ∆hmax(2ξij - 1) 1406
(5)
where ∆Ebs is the change in energy of the elastic bath. Using eq 7, ∆Ebs can be written as 1 ∆Ebs ) Cbijkl[V′′lk′ji - Vlkji] 2
(10)
where Cbijlk is the specified elasticity tensor of the elastic bath and ′ij is the strain tensor of the proposed move. The simulation is then carried out as outlined in section 2. The isothermal elastic constants of the composite material Nano Lett., Vol. 3, No. 10, 2003
Figure 1. Density as a function of temperature for the model polymer. The apparent glass transition temperature is Tg ) 0.5.
Figure 2. One component of the strain tensor, 33, vs MC step for the bulk amorphous polymer at T ) 0.8 with and without a positive elastic bath.
including the elastic bath, Ctijkl, can be calculated using eq 6. An apparent elasticity tensor for the material is then determined from Cijkl ) Ctijkl - Cbijkl
(11)
Note that the volume fraction of both the bath and the system of interest is unity. 4. Bulk Properties. In this section we apply the methods described in sections 2 and 3 to a model polymer system in the bulk. Each chain in our model polymer is made up of 16 interaction sites connected by Hookean springs. All sites interact via a 12-6 Lennard-Jones potential truncated at r ) 2.5σLJ and is given by uij ) 4LJ
[( ) ( ) ] σLJ rij
12
-
σLJ rij
6
(12)
where σLJ is the site diameter and LJ is the attractive energy between two sites at their equilibrium separation. Each interaction site is meant to represent a persistence length along the polymer backbone. In what follows, all values are reported in dimensionless Lennard-Jones units. This simple coarse-grained model is adopted in this initial study (over an explicit atomistic model) because it permits exploration of relevant time and length scales in reasonable amounts of computer time. Each simulation box contains 56 chains and the usual periodic boundary conditions are used. All simulations are conducted at an external pressure of 2.0. Note that at low temperatures, the density of the bulk system corresponds to the density in the interior of the nanoscopic structures discussed later in this work. To use the methods outlined in sections 2 and 3, a zerostress reference state must first be calculated. To do this, the polymer system is equilibrated in the NPT ensemble at each temperature in order to dissipate any residual, nonisotropic stresses and to determine the average density at each temperature. The reference box shape, 〈hij〉, for use in the NσT simulations, is then determined from the equilibrium Nano Lett., Vol. 3, No. 10, 2003
Figure 3. Contribution to the modulus of the composite system from the bulk amorphous polymer as a function of temperature. The glass transition temperature is at Tg ) 0.5. (see Figure 1)
density at each temperature. Figure 1 shows the density as a function of temperature for the model polymer considered here. An apparent glass transition temperature, Tg, can be associated with the temperature at which the slope of the density versus temperature plot changes. From Figure 1 we find Tg = 0.5. Simulations are then conducted in the NσT ensemble at an external pressure of 2.0 with zero applied stress, σij ) 0. Figure 2 shows one component of the strain tensor, 33, as a function of MC steps at T ) 0.8 with and without a positive elastic bath. The values µb ) 50 and λb ) 50 were used in eq 8. Figure 3 shows the contribution of the polymer to the modulus of the system as a function of temperature in the bulk. This modulus is given by E)
(3λ + 2µ)µ λ+µ
(13)
where λ and µ are the isotropic Lame´ coefficients of the polymer. The simulated, apparent modulus exhibits a marked increase with decreasing temperature, which can also be associated with an apparent glass transition temperature. The 1407
glass transition inferred from the mechanical properties (Tg ) 0.5) is consistent with that deduced from the thermal expansion curve (Figure 1). 5. Nanoscopic Polymer Structures. The mechanical collapse of densely packed photoresist structures has been attributed to strong capillary forces present during the drying of rinse liquids in wet chemical processing.17 Recent experimental efforts have been aimed at characterizing the mechanical stability of long photoresist lines.6,18 The nanoscopic structures investigated in this work consist of long rectangular polymer lines or plates standing on a substrate. These lines are infinitely long in one direction (through periodic boundary conditions), and finite in the other two directions. Three of the faces of the lines are “free” surfaces exposed to vacuum and the fourth is in contact with a substrate. The substrate consists of a face-centered-cubic (FCC) arrangement of interaction sites. Each monomer unit of the polymer interacts with each site of the substrate. Note, however, that throughout the simulation the positions of substrate sites are held constant. To generate initial configurations for the features of interest, thin polymer films were first formed using an expanded grand canonical ensemble (EGCE) technique on a lattice.19 Each simulation was started on an empty lattice. Chains were inserted into and deleted from the system at random. Other trial moves such as reptation, kink-jump, and crankshaft were performed to help equilibrate the system. To create periodic and nonperiodic boundaries selectively, hard surfaces were placed in the locations where periodicity was not desired. The equilibrium configurations prepared on the lattice were subsequently turned into a continuum (off-lattice) model by growing the lattice sites into monomer units having squarewell interactions. The films were then further equilibrated using discontinuous molecular dynamics simulations (DMD).8 DMD was chosen as the preferred method to equilibrate the thin films because of its computational efficiency. The nanoscopic structures (polymer lines) were then “cut” from the films by selectively deleting entire polymer chains from the system. This is analogous to using photolithography with a chemically amplified photoresist system. The manner in which the structure is formed from the film is likely to affect its final overall elastic behavior. A detailed description of the preparation procedure has been presented in the literature.8 The substrate was then introduced into the system and the resulting nanoscopic structures were equilibrated further on the substrate using DMD. The smallest structure considered in this work contains 576 monomer units, and the largest feature is made up of 28,224 monomer units (1764 chains). It is important to point out that extensive simulations of this model have shown that the glass transition of freestanding films is lower than that of the bulk material. For a film of thickness 9 σ (which is comparable to the smallest nanostructure considered in this work), the glass tramsition of the film is only about 7% below that of the bulk. All of the simulations presented in this work were conducted well below the glass transition (at ≈ 1/2 of Tg of the bulk polymer), thereby ensuring that the elastic properties reported 1408
Figure 4. Compression deformation simulations. A vertical force is applied downward on the top of the structure.
here are influenced by the proximity of a glass transition, even in ultrasmall structures. Finally, the interaction between monomer units was switched from square-well to the Lennard-Jones potential. Harmonic springs were inserted between adjacent monomer units in place of the flexible strings. The structures were then equilibrated further on the substrate using traditional isothermal molecular dynamics techniques. 5.1. Virtual Deformation Simulations. Explicit “deformation” simulations were implemented in a first attempt to determine the apparent modulus, E, of the nanoscopic polymeric structures prepared above. The rationale behind this approach is that they mimic experimental atomic-force microscopy measurements, and continuum-mechanics analyses which are generally used in conjunction with such experiments to deduce mechanical constants from deformation data. The deformation simulations were preformed using tradition isothermal MD. Although bending a structure from the side represents an actual deformation process present during the collapse of photoresist lines, a more direct method of calculating the modulus is to compress the feature from the top as shown in Figure 4. In this case the modulus is given by E)
σ
(14)
where the applied compressive stress is given by σ ) F/A and the resulting strain is ) ∆L/L0. Here F is the total applied force, and A is the cross sectional area of the structure; L0 is the unstrained height, and ∆L is the change in height. Figure 5 shows the applied stress, in units of LJ/σ3LJ as a function of strain for the compression deformation simulations. The modulus is taken as the slope of the linear regime. A minimum of four points along the stress-strain curve (Figure 5) in the linear regime was used to calculate the modulus. 5.2. Strain Fluctuations. The elastic moduli of the nanoscopic features were also determined using the strain fluctuation and elastic bath methods. These methods are described in detail in a separate publication20 and in sections 2 and 3. Nano Lett., Vol. 3, No. 10, 2003
Figure 5. Stress vs strain for three feature sizes under compression. The modulus is taken as the slope of the linear fits.
Figure 6. The modulus of the feature as a percentage of the bulk value as a function of line width determined from deformations (O) and strain fluctuations (b).
An elastic bath with positive elastic constants is used and the strain fluctuations are small. By calculating the entire elasticity tensor, the assumption that the features are isotropic need not be made. We do, however, assume that the features are homogeneous and report the elastic constants of the entire feature. It is noted here that these systems are likely to be inhomogeneous at “atomic” length scales, and a local definition of elastic constants may eventually provide some useful insights into the behavior of size-dependent mechanical properties. 5.3. Results. The aspect ratio (L0/w) was held constant at approximately 4 for all features. The temperature was set to T ) 0.2, well below the apparent glass transition temperature for this system Tg ∼ 0.5 in the bulk. Figure 6 shows the modulus as a percentage of the bulk value for the nanoscopic features determined from the compression deformation simulations. The modulus of the largest feature is nearly 100% of the bulk value. With decreasing feature size, the modulus decreases from the bulk value. It was not possible to simulate features with line widths larger than that shown due to CPU time requirements. Simulations of smaller features were difficult as well because these are mechanically unstable. Figure 6 also shows the contribution to the modulus from the polymer as a percentage of the bulk value for nanoscopic Nano Lett., Vol. 3, No. 10, 2003
Figure 7. The modulus as a percentage of the bulk value depending on the direction in which the measurement is taken. The periodic direction is perpendicular to the page.
structures determined from strain fluctuations. The same features from before were used here. The Lame´ coefficients of the elastic bath were set to λb ) µb ) 50. As before, the largest feature considered here has a modulus less than but close to the bulk value. The modulus decreases with decreasing line width. For the smallest feature studied here, the modulus is very nearly zero. The elastic bath technique permits study of both larger and smaller features than was possible with the explicit deformation method. Calculating the elastic constants from the strain fluctuations requires only a single simulation. Furthermore, the calculation converges faster than in the virtual deformation simulations. It is reassuring that both methods for estimating the apparent modulus of the features are in quantitative agreement. The virtual deformation technique uses equations from continuum mechanics (eq 14) and the assumption that the response to an applied stress is linear. The agreement of the two methods suggests that the features studied here can be assumed to follow a continuum-like behavior. However, the apparent size dependence of the modulus requires a molecular-level description not available in traditional continuum treatments. Both methods used in this study treat the structures as homogeneous materials. It is likely that these structures, being quenched glasses, have properties (e.g., density, elasticity, surface tension, etc.) which are position dependent. It has been shown that in thin crystalline films the shear modulus is lower near the surface than in the bulk.21,22 It is unknown wether only the regions near the free surfaces contribute to the overall size dependence or if the entire structure becomes mechanically weaker with decreasing size. By correctly accounting for the position dependent elastic constants it may be possible to use continuum mechanics theories to predict the mechanical behavior of nanoscopic structures.23,24 At this point the origin of the size dependence on the modulus is unclear. Figure 7 shows the modulus as a percentage of the bulk property value depending on the direction in which the measurement is taken. These correspond to three different components of the elasticity tensor: Cxxxx, Cyyyy, and Czzzz. 1409
An important finding of this work has been to show that the elastic properties of polymeric glassy nanoscopic structures become highly anisotropic. For the smallest feature considered here, the modulus in the infinite or periodic direction (i.e., perpendicular to the direction of the page) is comparable to that of the bulk material. The vertical direction (normal to the substrate) has an intermediate modulus, about onethird of the bulk value, and the side-normal direction has the smallest elastic constant, only about a tenth of the bulk value. The vertical-direction modulus is the same as that calculated in the compression-deformation simulation. This result is physically intuitive, in that the nanostructures lose their mechanical rigidity only in the directions subject to confinement (the width and the height). To the best of our knowledge, however, these calculations represent the first direct observation of this effect in polymeric nanostructures. 6. Conclusion. Literature simulations of deformation experiments have shown that the bending moduli of nanoscopic polymeric structures are significantly lower than those of the corresponding bulk material. Those simulations, however, relied on explicit deformations in which bendingplate equations from continuum mechanics were assumed to be valid at ultrasmall length scales. Furthermore, they were conducted on model polymers interacting through a discontinuous square-well potential. In this work a Lennard-Jones potential was used to describe interactions between polymer monomers. The use of continuum-mechanics formulas was circumvented by direct calculation of the elastic constants of nanoscopic features from fluctuations of the strain. It was found that deformation calculations and strain-fluctuation simulations give consistent results. The latter formalism is more satisfactory in that it does not invoke many of the assumptions implicit in the deformation approach. The strain-fluctuation simulations revealed that the elastic constants of polymeric nanostructures are no longer isotropic; below some length scale, they become dependent on the size and geometry of the sample. In the “smaller” or confined directions of nanostructures, an isothermal, size-induced plasticization was observed well below the glass transition temperature. As the critical dimensions of nanostructures decrease, the surface-to-volume ratio increases. Eventually, the amount of material at the surface becomes a significant portion of the entire structure and begins to effect its overall properties. The free surfaces apparently play a major role in the deviations at nanometer length scales from bulk behavior. The structures studied here have dimensions that are on the order of the characteristic size of the constituent molecules; the radius of gyration of the model polymer is about 3σLJ. The fact that the apparent elastic constants are anisotropic in nanoscopic structures but isotropic in the amorphous bulk indicates that the overall mechanical behavior of the structures depends on their shape or geometry. This suggests that the nanoscopic structures may be mechanically inhomogeneous. Because much of the photochem-
1410
istry used in nanolithography was originally optimized assuming a homogeneous continuum, these results also suggest that other aspects of manufacturing nanometer-sized devices may behave unexpectedly for current photoresist materials. Experiments and simulations have shown that, depending on the substrate, the glass transition temperature (Tg) of ultrathin films can be higher or lower than that of the bulk.5,7,25 That work suggests that by “tuning” the interaction with the substrate the size dependence of the apparent modulus may be controlled, thereby providing a mechanism to prevent the mechanical collapse of nanoscopic photoresist structures. We are currently pursuing simulations of nanoscopic structures on different susbtrates; the results will be presented in a future publication. Acknowledgment. The authors would like to thank the Semiconductor Research Corporation for financial support. References (1) Streck, C.; Melnichenko, Y. B.; Richert, R. Phys. ReV. B 1996, 53, 5341-5347. (2) Wendt, H.; Richert, R. J. Phys.: Condens. Matter 1999, 11, A199206. (3) Forrest, J. A.; Dalnoki-Veress, K.; Dutcher, J. R. Phys. ReV. E 1997, 56, 5705-5716. (4) SIA, The International Technology Road Map for Semicondutors”, http://public.itrs.net. (5) Fryer, D. S.; Peters, R. D.; Kim, E. J.; Tomaszewski, J. E.; de Pablo, J. J.; Nealey, P. F.; White, C. C.; Wu, W. L. Macromolecules 2001, 34, 5627-5634. (6) Cao, H. B.; Nealey, P. F.; Domke, W. D. J. Vac. Sci. Technol. B 2000, 18, 3303-3307. (7) Fryer, D. S.; de Pablo, J. J. J. Vac. Sci. Technol. 2000, 18, 33763380. (8) Boehme, T. R.; de Pablo, J. J. J. Chem. Phys. 2002, 116, 99399951. (9) Kaul, A. D.; Gangwal, A.; Wadhwa, S. S. Curr. Sci. 1999, 76, 15611566. (10) Kawai, A. J. Vac. Sci. Technol. B 1999, 17, 1090-1093. (11) Laudahn, U.; Fa¨hler, S.; Krebs, H. U.; Pundt, A.; Bicker, M.; v. Hu¨lsen, U.; Geyer, U. App. Phys. Lett. 1999, 74, 647-649. (12) Fay, P. J.; Ray, J. R. Phys. ReV. A 1992, 46, 4645-4649. (13) Allen, M. P.; Tildesley, D. J. Computer simulation of liquids; Oxford University Press Inc.: Oxford, New York, 1987. (14) Parrinello, M.; Rahman, A. J. Chem. Phys. 1982, 76, 2662-2666. (15) Ray, J. R.; Rahman, A. J. Chem. Phys. 1985, 82, 4243-4247. (16) Landau, L. D.; Lifshitz, E. M. Theory of Elasticity; Pergamon Press Ltd.: Oxford, New York, second ed.; 1970. (17) Tanaka, T.; Morigami, M.; Atoda, N. Jpn. J. Appl. Phys. 1993, 134, 16-30. (18) Stoykovich, M. P.; Cao, H. B.; Yoshimoto, K.; Ocola, L. E.; Nealey, P. F. AdV. Mater. 2003, 15, 1180-1184. (19) Escobedo, F. A.; de Pablo, J. J. J. Chem. Phys. 1996, 105, 43914394. (20) Van Workum, K.; de Pablo, J. J. Phys. ReV. E 2002, 67, 011505. (21) v. d. Eerden, J. P.; Knops, H. J. F.; Roos, A. J. Chem. Phys. 1992, 96, 714-720. (22) Van Workum, K.; de Pablo, J. J. Phys. ReV. E 2002, 67, 031601. (23) Miller, R. E.; Shenoy, V. B. Nanotechnology 2000, 11, 139-147. (24) Gurtin, M. E.; Murdoch, A. I. Arch. Ration. Mech. Anal. 1975, 57, 291-323. (25) Fryer, D. S.; Nealey, P. F.; de Pablo, J. J. Macromolecules 2000, 33, 6439-6447.
NL034458L
Nano Lett., Vol. 3, No. 10, 2003
