A Fully Coupled Computational Model of the Silylation Process

SANDIA REPORT SAND99-8225 Unlimited Release Prjnted February 1999 r F J A Fully C@pled’Computational Y4 th&~ilylatio-n Probess , // ,$/ --’c \ f”/ ~.\., y+..~y . .... 7 --- “+Lf- ;:~(= .Jy ‘ ,,, ‘- ... ‘\ -Y . -+ W. S. Winters~G.,H. Evans,.R. “S~”-LarsonfV. C. Prantil /’ .. /“” ‘x.. / .’ ,, -...+ ,,, “ /; . .. / Prepared by Sandia %onai Laboratories o /“ Y Albuquerque, New Mexico 87185 and Lwermore, California 94550 / / Sandia is a multiprogrsm l~boratory operated by Sandi{Corporation, a Lookheed Martin Company, for the United States Department of Energy under Contract’DE-AC(J4-94AL85000. ./ ,’ / 1 z~-d-for public releasq further dissemination unlimited. Laboratories / I / / Model of Issuedby SandiaNationalLaboratories,operatedfor the UnitedStatesDepartment of Energyby Sandia Corporation. NOTICE: This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government, nor any agency thereof, nor any of their employees, nor any of their contractors,subcontractors,or their employees,make any warranty,express or implied, or assume any legal liability or responsibility for the accuracy, completeness,or usefulness of any information, apparatus,product, or process disclosed,or represent that its use would not inhinge privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise, does not necessarily constituteor implyits endorsement,recommendation,or fiavoringby the United States Government, any agency thereof, or any of their contractors or subcontractors. The views and opinions expressed herein do not necessarily state or reflect those of the United States Government,any agency thereof, or any of their contractors. Printed in the United States of America. This report has been reproduced directlyhorn the best available copy. Avadableto DOE and DOE contractorshorn Officeof Scientificand TechnicalInformation P.O. BOX 62 Oak Ridge,TN 37831 Pricesavailablefrom (615) 576-8401,FTS 626-8401 Availableto the publictiom NationalTechnicalInformationService U.S. Departmentof Commerce 5285PortROydRd Springfield,VA 22161 NTISprice codes PrintedCOPY A06 Microfichecopy AO1 DISCLAIMER Portions of this document may be illegible’ in electronic image products. Images are produced from the best available original document. SA.ND99-8225 Unlimited Release Printed February 1999 A Fully Coupled Computational Model of the Silylation Process W. S. Winters, Computational G. H. Evans Reactive and R. S. Larson Processes Sandia National Laboratories Department / California V. C. Prantil Solid and Material Mechanics Department Sandia National Laboratories / California . This report documents the development of a new finite element model of the positive tone silylation process. Model development makes use of preexisting Sandia technology used to describe coupled thermal-mechanicalbehavior in deforming metals. Material properties and constitutive models were obtained horn the literature. The model is two-dimensional and transient and focuseson the part of the lithographyprocessin which crosslinkedand uncrosslinkedresist is exposed to a gaseous silylation agent. The model accounts for the combined effects of mass transport (diffusion of silylation agent and reaction product), chemical reaction resulting in the uptake of silicon and material swelling, the generation of stresses,and the resultingmaterialmotion. The influenceof stresson diEuSon and reaction rates is also included. Both Fickian and case II -Ion models have been incorporated. The model provides for the appropriate mass transport and momentumboundary conditions and couples the behavior (stress/strain) of uncrossli,nkedand crosslinkedmaterials as well including the underlying device topology. Finite element mesh generation, problem setup, and post processing of computed results is sufficiently mature to permit investigation of a broad parameter space which includes material properties and geometry issues (e.g., alternate crosslinlcingdistributions, resist thickness, etc.). The 2D transient model has been validated for simple geometries using independent materialpoint simulationsand onedimensional transient simulations. A complete 2D transient simulation of the silylation process is presented and discussed. Recommendations for future work are presented. 3/4 —. ._—.— 9 1 Nomenclature Greek symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.2 Subscripts and superscripts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 LI 13 2 Introduction 2.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2 Earlier Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.3 Report Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 17 3 Problem Description 3.1 TSI Processing Steps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.2 Silylation Modeling Considerations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 21 4 Model Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 4.1.1 Species Mass Conservation Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 4.1.2 Unconstrained Volume Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 4.1.3 Chemical Production Rates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 4.1.4 DiffuSon Flux Expressions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 4.1.5 Initial and Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 4.2 Momentum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 4.3 Polymer Material Response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 dChemistry . 4.1 Mass Tknsferan . . . . 4.3.1 PhysicalDeformationMechanisms ., . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 SegmentRotation . . . . . . . . . . . . . . . . . . . . 30 4.3.4 Evolution Law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.3.5 Resistance to Chzmgesin ConfigurationalEntropy . . . . . . . . . . . . . . . . . . . . 31 4.3.2 Kinematics. . . . . . . . . . . . . 4.3.3 Inelastic Flow Rule for Pol~er 33 5 Solution Method 5.1 Mass Transfer and Chemistry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 5.2 Momentum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 5.3 Pol~er 41 Material Response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 5.4 Tme Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 5.5 Mesh Generation and Postprocessing 46 5.3.1 Coupled Newton-Raphson Algorithm.. 5.3.2 Radial Retmn Algorith m....... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 6 Results 47 6.1 47 Constitutive Model Material Point Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1.1hiaxialExtension.... 6.1.2 Elastic Swehg. 6.2 Comparison oflDand2D 6.3 Transient2DSilylation -Notid Concluding A APPENDIX . . . . . . . . . . . . . . . . . . . . . . . . . . ...48 Pree Surface Swelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ...48 6.1.31dealized 7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ...49 Model Results . . . . . . . . . . . . . . . . . . . . . . . . . . ...49 Cdtiation Remarks - A One-Dimensional . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 57 Model 59 6 -.. .— . —---- -.—- . .. List of Figures 1 TSIprocessing steps. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 2 Cross-sectionalplasma stain of silylation prior to etch. . . . . . . . . . . . . . . . . . . . . . . 66 3 Arbitrary moving control volume. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4 A two-dimensionalquadrilateralcontrol volume centered about point P. . . . . . . . . . . . . 68 5 Control volurneand fluxes atresist/gas interface. . . . . . . . . . . . . . . . . . . . . . . . . . 69 6 The Eyring dashpot element (upper figure) and the Langevin spring element (lower figure). . 70 7 Dashpot response. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 8 Langevin spring response. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 9 Four node quadrilateralelementin global coordinate system. . . . . . . . . . . . . . . . . . . 73 10 Four node quadrilateralelementin local isopsmmetric coordinate system. . . . . . . . . . - . 74 11 A four node quadrilateralelementshowing sub control voluinesand integration points for the CVFEM. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 12 Response inunhxia lextension . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 13 Rate dependence of flow strength in Uniaxialextension. . . . . . . . . . . . . . . . . . . . . . 77 14 Softening response asafunction ofmodulus inunimdal extension. . - . . . . . . . . . . - . . 78 15 Stexlystate flowstrength responsein uniaxialextension. . . . . . . . . . . . . . . . . . . . . 79 16 Lzwgestrain unisxial extension simulation Comparison of material point simulationand 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 17 Response in freesurface swelling. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 18 Strain softening response in freesurface swelling. . . . . . . . . . . . . . . . . . . . . . . . . . 82 model 19 Freesurface expansionforidealized elastic response. . . . . . . . . . . . . . . . . . . . . . . . 83 20 Large strain swelling simulation Comparison of materialpoint simulation and 2D model. . . 84 21 Geometry forl-D solution comparisons. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 22 1-D and 2-D model predictions for layer thiclmessas a function of time. . . . . . . . . . . . . 86 23 D~tributions of species concentrationsand stress at three difEerenttimes for v = 0.5. . . . . . 87 24 Computational mesh for the nominal 2D silylation simulation.. . . . . . . . . . . . . . . . . . 88 25 2D material swellingwith superimposedconcentration proiilesof the silylation agent (species of time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 26 2D concentration profiles for species B, C, D, and E at 60 seconds. . . . . . . . . . . . . - . . 90 27 2Dstress distributions at60seconds. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 28 2Ddeformed meshat60seconds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 29 Concentration profiles for silylating agent P in one-dimensionalmodel for 6 = 0.01, /3 = 1, A)asafunction Tr=l, 30 ~=2, ~= l,v=0.5, aud~= Concentration profiles for unreztd ~r=l, a=2, ~=1, v=0.5, l. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 solid Q in one-dimensionalmodel for 6’ = 0.01, P = 1, and A=l. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 94 31 Concentration profiles for reacted but unexpanded solid U in one-dimensionalmodel for 9 = 0.01, #=1, 32 ~d~=l. . . . . . . . . . . . . . . . . . . . . . . . ~r=l, a=2, ~=1, z/=0.5, and~= l. . . . . . . . . . . . . . . . . . . . . . . . . . . l. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 97 l . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ...98 Progress of swellingin one-dimensional model for 6’= 0.01, @ = 1, ~r = 1, a = 2, Y = 1, and Ai l. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 95 Stress profiles in one-dimensionalmodel for O = 0.01, 6 = 1, ~, = 1, a = 2, -y = 1, v = 0.5, and~= 35 v=0.5, Velocity proilles in one-dimensionalmodel for O= 0.01, @ = 1, r, = 1, a = 2, ~ = 1, v = 0.5, and~= 34 cy=2, ~=1, Concentration profilesfor reacted zmdexpanded solid E in one-dimensionalmodel for r3= 0.01, ~=1, 33 ~r=l, Overall progress ofswelling inon+dixnensional modelfor O = 0.01, /3 = 1, Tr = 1, a = 2, v=l, ad~=l . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ...100 8 99 1 Nomenclature a A A A.(t) B B c c, c Q D D D Do E F Fi,j {F} f h H IJI ji Kg k k. k8 y] i m mip n N ~1 P“ Q q R RHS R s s t At u Ur UI %1 us Vi v Vi 61 mean molecular radius area,proportionality constant silylatingreagent surface area of arbitrary control volume (function of time) element geometry parameterdefined by Equation (115); amine product of silylatingreaction given by Equation (1) unexposed, unsilylatedpolymer resistmaterial rubbery modulus elastic modulus tensor molar concentration of species i (moles i / volume) unexpanded, silylated resistmaterial deformation rate tensor diEusioncoefficient (cf. Equations (21) and (22)) constantpart of diffusion coefficient (cf. Equations (21) and (22)) expanded, silylated resistmaterial body force vector, deformation gradient tensor nodal body force components where i = z, y and j = 1,2,3,4 column vector of body force components defied by Equation (86) density scaling factor hardening modulus Henry’s constant element Jacobian defined by Equation (77) mass *ion flux of species i mass transfer coefficient Boltman’s constant rate constant for polymer relaxation reaction (Equation (2)) silylation rate constant for silylation reaction (Equation (l)) mass matrix molecular weight of species i mass mass of species i in control volume centered about point p unit normal vector stressdirection elementbfiear shape function deilned by Equation (57) where 1=1,2,3,4 scalar material constitutive parameter rotation tensor scalar material constitutiveparameter yield limit right hand side rotation tensor flow resistance path along a surface time time step right stretch tensor x component of displacementat node 1,1 = 1,2,3,4 x component of velocity at node I, I = 1,2,3,4 x component of acceleration at node 1,1= 1,2,3,4 y component of displacementat node I, 1=1,2,3,4 velocity of species i left stretch tensor WfuAon velocity of species i y component of velocity at node I, I = 1,2,3,4 9 VI v . ;} v.(t) Vp(t) w w w~b w z x z~ Y ~ YI Y 1.1 a q [A] @ y component of acceleration at node I, I = 1,2,3,4 velocity vector (mass average) accelerationvector column vector of acceleration components defined by Equation (85) volume volume of arbitrary control volume (function of time) volume of control volume centered about point P (fnnction of time) spin tensor finite element weightingfunction arb@iry velocity of control volume coefficient in exponential *ion expressions (cj Equations (21) and (22)) Cartesiancoordinate x position vector x coordinate of node I, I = 1,2,3,4 cartesiancoordinate y mass fraction of species i y coordinate of node 1, I = 1,2,3,4 yield function Greek symbols coefficient of thermal expansion volume fraction of species i coefficient matrix defined by Equation (88) nodal dependent variable defined by Equation (54), silylation swelling density intrinsicdensity of species i (mass i / volume i), stresstensor column vector of stress components defined by Equation (87) x component of the traction vector on a plane whose normal points in the x direction y component of the traction vector on a plaae whose normal points in the y direction z component of the traction vector on a plane whose normal points in the z direction x component of the traction vector on a plane whose normal points in the y direction (aZV= cry=) local isoparametnc coordinate local isopamnetnc coordinate, angle of molecular chain rotation mass production rate of species i plastic strain rate deviatoric stresstensor magnitude absolute temperature Lard elastic constants Poisson’s ratio Langevin back stresstensor, rate dependent yield stress — 1.2 m d ._. . . . ... Subscripts and superscripts momentumquantity filon quantity point P at the control volume center 9,s evaluatedin the gas at the surface g,co evaluatedin the gas far from the interface species index i point P at the control volume center P r,a evaluatedin the polymer resistat the surface elastic e plsstic P P 11/12 2 Introduction This report describes a new two-dimensional transient computational model of a positive tone silylation process. The model accounts for multispecies dif?kion of the silylation and product gases, the silylation chemical reaction, the subsequent swelling and movement of the resist and adjacent materials, and the evolution of the stress state including the kinematics for finite deformations. Although the model is tw~ dimensional, its extension to three dimensions is straightforward. Material properties and constants were estimatedfrom the literature. The abtity of the model to predict accurate quantitativebehavior will depend on the collection of more precise material property data. 2.1 Background Photolithography and its related chemical processing steps are ~ected to have an essential role in the development of next-generation semiconductor devices. Extreme UltravioletLithography (EUVL) has been demonstrated as a viable candidate for fabrication of integrated circuits having feature sizes of 130nm or less. Due to the strongly absorbing nature of 13.4nmradiation, thin layer imaging (TLI) resist technology is utilized by EUVL. In TLI, an imaging mask layer is placed over a bottom layer photoresist. The top layer is then irradiated causing the photochemical transfer of the mask image to the top surface of the photoresist. Henderson et cd. [1] discuss severalTLI processes. The imaged pattern in the top thin layer of resist is then pattern transferee into the remaining resist thicknessusing an anisotropic plasma etch step. The post-exposure silylation process, also referred to as top surface imaging (TSI), is one of several viable processesbeing proposed for TLI (see e.g. reference [1]). Like other TLI methods, TSI is an excellent match to EUVL because it allows thick resist layers to be patterned despite the limited penetration depth of the incident radiation (see e.g. reference[l]). The critical steps of TSI include exposure, post-exposure bake, silylation, and etch development. These steps have an influence on process sensitivity, contrast, and resolution. Depending on the chemistry of the photoresist, either the exposed(negativetone) or unexposed(positivetone) regionscan be selectively silylated. During the silylation step, gaseous ami,nosilanesare absorbed at the resistsurface and dif7useinto the imaging layer. A reaction takes place in the layer that resultsin an uptake of silicon, localized increases in volume, and the evolution of a product gas. Volume increases on the order of one hundred percent are typical which results in the development of complex stress states in the resist material. These volume increasescan influencethe effectivenessof subsequentprmessing steps (e.g. etch development). Attempts to control swellingof the resist during silylation are discmsed by Han et al. [2]. Researcherswho have examined the behavior of materialsundergoing TSI include La Tulipe et d. [3], Horn, et uZ.[4], Itani, et al. [5], Glezos et al. [6], Whelan, et uZ.[?, Zuniga and Neureuther [8], Ocola and Cerrina [9], and Henderson et UL[1]. 13 2.2 Earlier Work Early silylation modeling efforts represented silylation solely as a diffusion/reaction process. Weiss and Goethals [10] utilized a modified version of the Deal-Grove oxidation model [11] to simulate the onedimensionalpropagation of a silylation front as a function of time for negative tone resists. Bauch et al. [12] studied the silylationprocess experimentallyand numericallyusinga two-dimensional-Ion model. They used an equilibrium balance between sorption, resorption and diffusion as a boundary condition at the gin/resist bo~dary. The influence of resist glass transition temperature on the diffusion mechanism was examined by Hartney et d. [13] ad Paniez et al. [14]. They observed that, while Fickian &fFusionmay be appropriate below the glass transition temperature, Case II -Ion is more likely above it. This point is discussed in more detail in Section 4. A comprehensive difl!usion-basedsilylation model was proposed by Pierrat [15]. His model accounted for the reaction of the silylating agent with the resist polymer and the relaxation rate of the polymer after reaction. The diffusion coefficient of the silylating agent was calculated as a function of the relaxed polymer chain concentration. The models described above do not account for swelling of the resist duringsilylation. Such swelling can increase diffusion path lengths, induce stresseswhich affect the mode of &ion (Fkkian versus Case 11), and result in shape changeswhich can influence subsequent processing steps. A more comprehensivemodel of the silylation process must couple the effects of diffusion and large deformations in the material. Winters and Mason [16, 17] proposed a computational.model for heat conduction (diffusion of heat) and material motion for the purpose of modeling coupled thermal-mechanicalphenomena. The silylation model proposed here is an extension of their technique. Dimitrienko [18] developed a filly coupled model to determine the effect of tilte deformationson internal heat and mass transfer in elastomer ablating materials. TO the authors’ knowledge, the only silylation model to date which couples -ion of the silylation agent to the subsequentgrowth of the silylated medium was proposed by Zuniga and Neureuther[19]. Their model accounted for the relaxation of the polymer during silylation, the increase in diffusion rate after silylation and a slowing of the local silylation reaction due to the stress induced by swelling. They chose a variable diffusion coefficient similar to Pierrat’s [15] to account for an increase of diffusion through the ailylated medium. A silylation reaction rate was proposed that depended on the change of state achieved in the resist upon exposure and was proportional to the number of bonding sites available for the silylating agent. The reaction rate was also influenced by the stress generated during silylation. A detailed description of the constitutive equations and kinematicsof the model was not given. 2.3 Report Overview In the sectionswhich follow, the detailsof a new two-dimensionaltransientsilylationmodel will be presented. Model formulation, numerical implementation issues, results, and recommendations for future work are 14 discussed. Section 3 &scribes the silylationprocess and its relationshipto other TSI processing steps. Formulation of the silylation model is presentedin Section 4 including all appropriate equations and boundary conditions. Solution of the coupled set of governingequations is discussedin Section 5. Section 6 presentsseveralmodel verification studies md resultsfrom a 2D silylation analysis. Concluding remarksand recommendations for future work are presented in Section 7. A one-dimensionaltransient silylation model based on the model of Zuniga and Neureuther [19] is documented in Appendix A. 15/16 .. -.. — 3 .-+. / : .’>s ...7,....., , ...,.. .. .fe _ . ... ., ~,,, . : . .:, . ....zK.L_ ...,... ___ A _ Problem Description This section will describe the TSI process with an emphasison the silylation step. The nomenclature and terminology introduced here will be used in the discussionpresentedin later sections. Much of what follows is a summmy of information presentedby Wheeler, et.ul. [20]. 3.1 TSI Processing Steps A schematic showing the TSI processing steps and structure cross-sectionsare shown in Figure 1. Figure 1 (a) shows the cross-section prior to the first processing step. The cross-section consists of three layers, the substrate (black), an intermediateprocessing layer (grey) and the top surface imaging layer (white). The processing layer, a hard baked resist typically 400-500 nrn thick, functions as a pharizing layer and antireflectivecoating for the device substrate. The imaging layer, a“photoactive resist, is typically 300-400 nm in thickness. Figure 1 (b) showsthe first processingstep which consistsof applying a mask and irradiating portions of the imaging layer with 13.4 nm radiation. This results in photoinduced acid generation beneath the exposed surface. The cross-section subjected to acid generation is depicted here as the area enclosed by a dashed line. In this scenario, mid generation does not occur through the entire thicknessof the imaging layer. The next step, Figure 1 (c), is the Post Exposure Bake (PEB). A typical PEB occurs for two minutes at a temperature of 44)0K. This heating process causes the acid bearing portions of the imaging layer to be highly crosslinkedand nearly impermeableto silylation agents. The figure implies uniform crosslinkingover the acid bearing region. The actual uniformity and coverage of cmsslinkedregion is uncertain. Figure 1 (d) shows the silylation step for the positive tone process. Silylation increasesthe silicon content in certain areas thereby increasingetch selectivity for subsequentprocessing steps. The top surface of the imaging layer is exposed to a gaseousaminodiaiiane(silylationagent). Ty@al processing conditions call for low pressures(20-100 torr) and intermediatetemperatures (300-500K) for approximately one minute. The silylating agent diEusespreferentiallythrough the unexposed regions of the imaging layer where a localized chemical reaction occurs causing an uptake of silicon, a change in volume, and the liberation of a product gas. The magnitude and shape of the volume expansion depends on the reaction/-ion process and the resultingstressstate in the material. Jn Figure 1 (d) the resultingailylationregion is shown in red. Shallow silylation depths of approximately 20 nm are generally observed in the exposed areas. Such unwanted silylation is believed to occur because of acid loss during PEB by volatilizationfrom the surface of the resist. A Chlorine/argonplasma “descum” process is used to remove the thin silylated layer in exposed areas. The 17 descurnstep removes all residuesand improvesthe reproducibfity and linearityof the process [20]. Figure 1 (e) shows the device cross-section after completion of the descum step. The final step in TSI is the oxygen plasma etch, shown here in Figure 1 (f). The presence of silicon in the unexposed portions of the imaging layer restrictsetch penetrations to the exposed areas. The objective in TSI is to transfer the irradiated mask image to a thin imaging layer so that subsequent processescan be restrictedto cross-sectionswhich lie beneath either the exposed or unexposed areas. Fignre 1 represents an ‘idealization” of the TSI process in that all processing steps result in process boundaries which are perfectly orthogonal to a “flat” top surface. In reaMy, the process boundaries may be far from orthogonal due to aerialimage aberrationssuch as light scattering, nonuniformacid generation, nonuniform cross-linking,multidimensional-Ion of the silylatingagent, multidimensionalgrowth and deformationof the layer materialsand multidimensionaletching. Figure 2 illustrat~ this point by showinga cross-sectional plasma stain made after the silylation process and prior to the plasma etch. The series of equally spaced “bumps” along the top portion of the silylated resistcorrespond to the unexposed portions where silylation and volume expansion have taken place. Note also the presence of an unwanted Siiylation “scum” layer in the exposed regions. 3.2 Silylation Modeling Considerations The silylationmodel documented here is intended to representthe part of the TSI process which occurs after the PEB and before the descum step, see e.g. Figure 1 (d). This model does not predict the location and distribution of the crosslinkingprofile prior to silylation. Rather, this profile representsan initial condition for silylation modeling computations. The model carIaccommodate any crosslinkingprofile provided it is known a priori. The precise location of the crosdinking boundary, the uniformity of crosslinking,and the effectivenessof the crosslinkedregion in inhibiting diffusion and reaction of the silylatingagent axenot well known. Zuniga and Neureuther[19]have postulated that the potential for silylationcan be calculatedas a function of polymer characteristics,processing and imaging conditions. In their modeling they assumedthat crosslirddng could be related to energy deposition during imaging. They utilized a numericalsimulation(see e.g. [21]) to compute the energy deposition and ultimately the crosslinkingprofile. While the present model is capable of treating =Ion mechanisms ranging from pure Fickian to more complex mechanismswhich depend on local stress Statesand reaction rates, the precise nature of diffusion during silylationis not welllmown. Furthermore,mod~ parametersSU&as *~on coefficientsand reaction rates have not been measured for most of the resistsunder consideration for TSI. The present silylation model couples the complex ei%cts of di.fhsion, volume change and finite deformation stress development for the silylation and crosslinkedportions of the imaging layer. Furthermore,the model includes the compliance of the planari.zinglayer and device substrate when necessary,to compute stresses and deformations over the entire cross-section of the device. Accurate model predictions require knowledge of mechanical material parametersthroughout the device cross-section. At the presenttime, many of these parameterscm only be estimated horn similar materials. Until precise experiments can be designed to determine diffusion and reaction mechanismsand to measure transport and material parameters, predictions from this or any other silylation model must be regarded as qualitative. Nevertheless,the present model should be extremely valuable in assessinga wide rage of silylation processing scenarios. 19/20 _.=. 4 ., .,., ~i,<f-k——— . , ,---,......... .-,-., --- -“2—— . .. . . Model Formulation The present silylation model couples the effects of silylating agent and reaction product mass transfer, silylation reaction, volume change, stressgeneration, and iinite deformation of the device cross-section. This section describes the formulation of the mass transport equations, the equation of motion (momentum), and the expression for unconstrained(stress-free)volumetricexpansion due to silylationchemistry. Also included is a discussion of boundary conditions and sub-models for materialconstitutive behavior and mass diflision. The formulation presented here is geometry independent. 4.1 Mass ‘lhnsfer and Chemistry The silylation process studied incorporatessilicon into the unexposed (and uncrosslinked)areas of the resist (typically a crescdnovolac) thereby renderingthose areasresistantto “subsequentetchin~ phenomena include diRMon, chemical reaction, and large (100%) volume changeof material. The process takesplace in a closed, isothermal chamber that is maintained at a temperature of 70-80 C and filled with a silylating gas (e.g., dimethylaminopentamethyldisiie, DMAPMDS) to a pressureof approximately 30 Torr (cf. reference [20]). In the model the silylating gas is assumed to be adsorbed onto the unexposed surface of the resist; once absorbed the silylating reagent (denoted by A) difl?usesinto the unexposed and unsilylated polymer resist material (denoted by C), reacts with it forming unexpanded, silylated resist (denoted by D) and releasing amines (denoted by l?) which diffnse through the material and are desorbed from the resist surhce. The silylating reaction considered in this study is given by A+ C%?+D (1) A second reaction (polymer relaxation), which accounts for the swellingof the resistdue to the mass addition included in the silylating reaction (1), is given by D%E (2) This chemical model has been used by Pierrat [15] and Zuniga and Neureuther [19]. from a material modeling point of view, it is desirable to representthe process in a Lagrangian reference frame, i.e., a reference frame moving with the deforming material. The derivation below accomplishes th@ however, mlon to take phe of the silylatingreagent,~ and the product species, 13,of the silylatingreaction is allowed across the system (Co-mg of species C and, afteI chemical reactions (1) and (2), species D and ~ as Well) boundaries. 21 4.1.1 Species Mass Conservation Equations The approach adopted here for species mass conservation is the arbitrary moving control volume approach of Whitaker [’2’2]where the control vohune velocity is W=b (cf Figure 3). This approach yields the following equation for mass conservation of species i : d ~ / v=(t) pYid-v + (&w = * (t) Pi(vi – w~b) . ndd – /Va(t) [ a o (3) When the velocity of species i equals the velocity of the control volume (vi = W=b) this equation becomes: d ~ I v.(t) pYidv – /Va(t) lilidv = o (4) Site we choose a condition that restricts the unsilylated and unexposed material C to remain within the moving volume, the control volume velocity, w, is not arbitrary; indeed it must equaJthe velocity of species C, vc. Reactions (1) and (2) form D and E, and these materials are also contained within the moving volume so that: (5) v~ =v~=v~=w. The velocity, w, of the control volume is in generalnot the mass averagevelocity, v. We deiine dfilon in the usual way relative to the mass average velocity, v: jj = pYi(vi – v) = pYivi (6) where v~=v~–v. is the W?u.sionvelocity of species i. Note that the sum of the mass -Ion E ji=o. Equations (5), (6), and (8) can be combined to yield: 22 (7) fluxes must be zero, that is: (8) v—w= j. +jB (9) Pc+PD+PE” The mass conservation equations of materialsC, D, and Z which remain within the moving control volume VP(t) (shown in Figure 4 as a two-dimensional quadrilateralcontrol volume centered about point P) are given (for mass production rates that me constant over the control volume) by: dmc, — = (#l&v-(t) dt (lo) dm. ~ dt (11) = tiDpvp(t) dmE ~ = (&,vp(t). (12) dt Note that although there can be -Ion of materials C, D, and E across the moving control volume boundary (e.g., V, s v. – v = w – v # O), there is au equal and opposite flux of these materialsconvected across the moving control volume boundary by the mass average velocity, v, resultingin zero net transport of materialsC, D, and E across the boundary. There is a flux of materialsA and B across the boundary of the moving control volumq the mass conservation equations for A and 1? are given by din., — dt = ti.pv~(t) – dmBP — dt = tiBpvp(q – /[ Ap(t) 1[Ap(t) pAf5A(jA+j.) j. + (Pcec + %% PB~i9(jA j. + (Pcec + Z&) +&) + PIXD + E&E) 1 (13) 1 (14) . ndA . ndA. The quantity j3iis the intriniic density of species i (gramsi /(cm3 i )); e; is the volume fraction of species i (volume i /volume mixture), and jA and j~ are the mass -Ion 4.1.2 Unconstrained fluxes of species A and l?, respectivdy. Volume Equation There can be significant volume change of the resistduring silylation. Estimatesmade by Wheeler et al. [20] of volume changes of a cre.ol novolac resist silylated completely using various ailylatingagents ranged horn 12% to 124%; these estimates were based on the chzmgein the number of atoms in the repeat unit of the polymer. 23 If the intrimic densities, ~i (grams i/cm3 i ), of the species are assumedto be constant and if the volume fractions of the species sum to one (~ ei = 1), then the species mass conservation Equations (10)-(14) can be added to give an explicit equation for the time rate of change of the volume of the moving control volume: +jB) 1 - ndd. (15) This equation describes the unconstrained volumetric expansion due to silylation. . 4.1.3 Chemical Production Rates The silylation reaction (1) assumesthat one molecule of the silylating agent A reacts with one repeat unit of the unexposed and uncrosslinkedpolymer C to form one molecule of product species 13 and one repeat unit of silylated, unexpanded polymer D. The relaxationreaction (2) impliesthere is a time scale associated with swellingof the silylated polymer ( [15], [19]). The mass production rates of the species are given by: k.p.pc 4A = –— M= _ k8pApcMB “ – M.MC (16) (17) (18) ~~ = k8pApcM. – kep= M.MC WE= kep.M. M. (19) (20) where the units of k~ are cm3/mole-see and the units of ke are see-l. Zuniga and Neureuther [19] allowed the silylation rate constant to depend exponentially on the stressin the resist. They referred to the work of Hm-tney[23]who postulated that a silylationfront propagatesinto the resistby a Case II diffusion mechanism ([24], [25]) in which the silylation rate depends exponentially on stresswhen that stress exceeds a critical swellingstress ([26]). We note that in the resultspresentedlater in this study the silylation rate constant k~ ifwarrantedthe presentmodel can be mod.i&d easily to allow a stress dependent is assumedto be constant; rate constant (cf. Appendix A). 24 4.1.4 DiEu&on Flux Expressions The mass diffusion fluxes discussedabove are given for species A and B by j. = ‘@v(?n./m) = ‘p~.e (wpE/ME)v(mA/m) (21) = –pDoe (w%f~’)V(m~/m) (22) and jB = –pDV(mB/m) Equation (21) is sixnilarto expressions given in [15] and [19]; for w = O Equations (21) and (22) are expressionsof Fick’s law (where for lack of data we have assumed e@al diffusion coefficientsfor .4 and B). Pierrat [15]proposed the exponential dependence of the diffuion coefficient on the expanded, silylated resist concentration c~, arguing that the filon coefficient should be larger in the silylated areas due to the swellingand decreasesin density and glass transitiontemperatureof the polymer that occur in the silylation process. Note that Crank [271assumed an exponential dependence of the diffurion coefficient on penetrant concentration to describe non-Fickian dif?usionin polymers. 4.1.5 Initial and Boundary Conditions The dependent variablesin Equations (10)-(15) me the massesof the species and the volume of the control volume; thus the initial and boundary conditions depend on volume. For example, for a control volume centered about point P (cf. Figure 4), miP(t) = pi=(t)VP(t). Since initial and boundary conditions are usually sptied in terms of concentration, the discussionhere will refer to (mass) concentratiorqconversion to species massestakes place in the computer code after discretizationwhen the volumes are available. The initial conditions for the species mass conservationequations in the region where si.lylationis allowed to occur state that only C is present: pc (x, O) is specified and p. (x, O) = p~(x, O) = p~(x, O) = P=(x, O) = O. Boundary conditions ~e required for species A and B. For all boundaries other thaa the gas/polymer interface, the boundary conditions ~e zero normal fluxes of species A and B: j~ . n = .L “ n = OAt the gas/polyner interface two types of boundary conditions are implemented. For the first type, a fhed concentration of A and zero concentration of 13are specifie& pA,,~ is specifiedand p~,,~ = O. For the second type, a steady-state mass balance that includes an equilibrium assumption at the surhe is the specifkd boundary condition for A; zero concentration of B is againspeciiied. This second type of boundary condition 25 for A is a flux boundary condition that can be explained with respect to the gas/polymer interface shown in Figure 5. As the thiclmess (A~) of the control volume containing the surface shown in Fignre 5 sbrinks to zero, the mass balance for the silylating agent A is: where hAr is the mass flux of silylating agent A at the surface into the resist and ?hA9 is the mass flux of silylating gas A into the surface from the gas, which can be given by: rnAg = Kg bAg,w – PA9,8)” (24) where Kg is the mass transfer coefficient (ii appropriate units), pAg,m is the partial pressureof A far from the interface, and pAg,~is the partial pressure of A in the gas phase at the gas/polymer interface. At the interface, equilibrium is assumed amd the partial pressurep~~,~ is related to the mass fraction YAr* (YA s p./p) of species A in the resist at the surface by (Middleman [28]): (25) where H is Henry’s constant (a function of temperature). Combining (24) and (25) yields: ( ?f2Ag= Kg p.g,w – ~KT,8 Mmixt.rer,8 MA ) (26) which gives an equation for the unknown mass fraction (or concentration) of A at the surface of the resist in terms of the parametersKg and H and the specified chamber pressurep~g,m. 4.2 Momentum The material motion is obtained iiom the conservationof momentum, i.e., a(p) —= t% V.C+F (27) where p, v, ISand F are the material density, velocity, stress, and body forces respectively. The body force is included here for completeness. The only applicable body force for the sily~lon problem is that due 26 ,.., -—. . .... . . __=zw.> ....< - .---AA , -A—— — —. —. to ~avity which is neglected in the numerical implimentation. The inertia term on the left hand side of Equation (27) is also negligiblesince the materialmotion is quasi-static. For reasonsexplained in Section 5, a portion of this term will be retainedto facilitate an explicit time integrationfor v, the principal dependent variable in Equation (27). The material velocity can be related to the stress through the deformation rate tensor. This relationship varies depending on the kinematics snd constitutive model employed (see Section 4.3). Boundary conditions imposed on the equation of motion are straightforward. At the free surface of the imaging layer, the velocity field is computed in such a way as to insure the following stress free boundary condition: (28) were n is the unit normal to the surihce. The interfacesbetweenthe imaginglayer, processinglayer, and devicesubstrate are assumedto be bonded together and hence part of the interior computed solution, i.e., there is no need to provide boundary conditions for material interfaces. At some depth sufEcientlyfar from the free surface,the materialis immobile suchthat the following boundary condition holds: V=o (29) Boundary conditions placed on the in-plane and lateral directions of the device cross-section depend on whether the analysisis three-dimensional,two-dimensional,or axisymmetnc. For the two-dimensionalplanar model developed here, symmetry boundary conditions are assumed. 4.3 Polymer Material Response This subsectiondiscussesthe polymer constitutive model. The physical deformation mechanisms,model kinematics, inelastic flow rule, evolution law, and backstressmodels are formulated in the sections which follow. 4.3.1 Physical Deformation Mechanisms To deform inelasticaUy,polymers must overcome two distinct internalresistances. Initially,the flow strength is controlled by an internal resistanceto deformation-induced rotation of the individual molecular chaius 27 that comprise the microstructure. This resistance corresponds to breaking bonds connecting randomly oriented nearest neighbor polymer chains. Subsequently,Aains undergo aiiinerotation and realign at larger strains. AS the polymer deforms into this more orderly state, an additional resistmce must now be overcome corresponding to this increasein conf@rationaJ entropy. Often, polymer resistsare described by a viscoelasticconstitutive law. Further,the recoverable elasticstrains are almost always small. As such, the material responseis well-representedby a viscous element in parallel with an elastic spring. Below the glass transition temperature,however, the viscosity of the inelasticelement is very large with a correspondinglylarge time constantfor relaxation of any of the inelasticstrainsl. So while a viscoelastic constitutive law is often used, the appment permanence of strain below the glass transition temperature makes a simpler viscoplastic law sufficient. 4.3.2 Kinematics Glassy polymers can, in general, exhibit both elasticand inelastic deformation. In what follows, we exercise the inelastic glassy polymer model developed by Boyce, Parks & Argon and described in detail in [30]. The formulation used here is restricted to small elastic strains while allowing for arbitrarily large inelastic deformation. The motion is prescribed by a multiplicativedecomposition of the total deformation gradient into elastic and plastic parts F = F“FP = V“RUP (30) where V and U are the left and right stretch tensors, respectively. Furtherrestricting the elastic deformation gradient to be symmetric F ●T =F”=V* (31) k.u’ (32) the plastic deformation gradient F’ = then describes the relaxed ccdguration upon elastic unloading to a stress free state without rotation. It follows that the velocity gradient may be written 1“Because the internal viscosity is relatively high, only a smaUpart of the deformation is recovered at normal temperaturesso that the plastic deformations have the appearance of permanence.” [29] 28 +V” D+ ’ wp [ 1V*-l L = @/”-l (33) where Dp and W“ are, respectively,the symmetric and skew parts of the plastic velocity gradient, known respectively as the plastic rate of deformation and the plastic spin. The Cauchy stress, U, is uniquely defined by the natural logarithm of the elastic deformation gradient 1 C“ : @{v”} c = det(V”) – (cwA@ + A@)I’j (34) where the volumetric terms are the thermal and silylation swelling , respectively, and the elastic modulus tensor is given by c“=2pII+AIt31 (35) where p and A are the elastic Lam6 constants and II and I are, respectively, the fourth and second order identity tensors. The swelling A@ is the net unconstrained volumetric expansion calculated from Equation (15) For small elastic strains,the development by Hoger [31] can be used to show that ~(V”)]” N D“ (36) where D“ is the elastic deformation rate tensor. Using the finite deformation kinematics outlined in Bammann & Johnson [32] and Bammann & Aifantis [33], a rate formulation for the Cauchy stress cam be written as &=u —Woe + aW” N C- : [lnV-]” = C- : D- (37) where D“=D– Dp–{aO+@}I (38) and w“ =W–wp 29 (39) The quantity 6 is the unconstrained silylation strain rate and is determinedfrom. (40) where ~ ia determined from Equation (15). For small elastic strains and, owing to the symmetric choice for F“ [30], (41) - zmd2 U= C”: D- 4.3.3 (42) Inelastic Flow Rule for Polymer Segment Rotation The initial resistance to polymer deformation has its origin in the restriction imposed on molecular chain rotation horn neighboring chains as depicted in Figure 6. Inelastic deformation commences once a free _ energy barrier to molecular mobility is surpassed. This barrier is overcome by thermally activated rotation of chain segments under stress. The result is a Boltzmann expression for the inelastic deformation rate magnitude P=Joexp[*{l-(~):}l (43) ‘r where ~0 is a reference strain rate, SOis the athermal flow resistance, i.e. the resistance at absolute zero, stressmagnitude, p and q are material psrmneters and A is a 0 is the absolute temperature, ~ is the 10CSJ proportiomdity constant given by ~ = 397r<2a3 16k (44) Here & is the net angle of molecular chain rotation, a is the mean molecular radius and k is Boltzrnann’s constant. 2NOTE: The approximation resolved in (41) indicates that W* w O. Under an arbitrary superposd rotation, Q(t), We would be given by We = Q(t)QT(t) thus guaranteeingthat the rate expression in (37) remains objective. 30 - A typiczdillustrationof this strain rate dependence on stressis given in F@re 7. In this example, an Eying dashpot element would exhibit negligible strain for stressesbelow the temperature-dependent yield, here approximately 25 – 30 x 107dynes/cm 2. Above such stresses,the molecular chain segmentrotation becomes progressivelyeasierresulting in a strain rate that increasesrapidly with stress. We assume an associative flow rule with the plastic deformation rate and deviatoric stress coaxial D’ = ~PN (45) (46) 4.3.4 Evolution Law The response of glassy polymers is characterized by stress relaxation after yield. The stress relaxation is modeled by allowingthe internalresistanceto decreasewith continuedstraining. The Vote-type phenomen~ logical softening evolution described in [30] is used here (47) where s is the current athermal deformation resistancein the Eyring dashpot, S*8is the steady state value at large strain and h is the softening modulus, i.e. the slope of the yield drop with respect to inelastic strain. The initial value of this athermalresistanceis given by o.077p —— ‘0– l–v (48) The differencebetweenthe initial yield and the saturation valueof the dashpot resistanceis the value of the total stress relaxation in uniaxid extension. In this way, the chamxt~lc time constant for the inelastic response is determined by the model parameters h and sS*. A @pie.al strain-induced yield drop during volumetric swellingwill be illustratedin the section describingindividual response modes. 4.3.5 Resistance to Changes in Configurational Entropy At very large strain, there is a second contribution to resistanceto inelastic deformation. Associated with large rotations of long chain moleculesis a change in configurationalentropy that corresponds to resistance 31 of the bulk polymer to large scale alignmentand ordering of the molecular chains. In [30], a highly nonlinear Langevin spring elementrepresentsthe locking behavior consistentwith a non-Gaussianstatisticalmechanics description of rubber elasticity. In this formulation, an internal back stress develops in the polymer chain network as illustrated in Figure 6. This internal stress represents the ever increasing stress locked in the iiigned structure, i.e. the stress beyond which the polymer must be loaded to further strain and align the chain network so that now the constitutive law becomes (a - /3)” =C”:D” (49) The back stress is assumedto remain coaxial with the left plastic stretch tensor, VP, with principaJvalues given by (50) where F’ = RUP = VPR (51) and IV is the number of rigid chain links between entanglements. C’R is the rubbery modulus, and L is the Langevin function given by (53) This locking behavior is shown in Figure 8 where the nonlinear behavior of the Langevin function is illustrated. Effectively, the alignment of polymer chfi requires increasingly larger stresses. Because the strength of the back stress is proportional to the rigid chain link density, this spring element provides an appealing means of modeling crosslinkedpolymers through its observed locking behavior. Because the swellingstrainshave been shown to be of order unity, we describe this element of the complete model for illustration. However, as the locking behavior may not be realized by silylating polymer resists and because it has no counterpartin wmehstic . materialmodels, this elementis not exercisedin simulations reported here. 32 5 Solution Method This section describes the numerical implementationof the silylationmodel for twcdirnensiona.1(2D) transient problems. The equationsformulated in Section 4 me reduced to a set of fist order ordinary difkrential equationa that can be explicitly integrated in time. The computational space is discretized into 2D Lagrangkm finite elementsand all conservation equations are solved simultaneouslyon the same moving mesh. The model has been formulated so that each element contains a volume of solid mass that may or may not undergo a chemical transformationdepending on the materialtype and whether it is crosslinked. During a simulation, no solid mass is permitted to cross the boundaries of an element. The element may, however, deform and grow in volume due to silylation. All mass transport across the boundaries of an element is due to difli.Monand convection of the silylation agent and rea,ctionproduct. Figure 9 shows a typical four node quadrilateralelementin the globid 2D planar (z – y) coordinate system. Each elementis assumedto have a uniform and constantthicknessof unity in the z direction (the plane strain assumption). The goal in finite element discretization is to obtain a set of ordinary difkential equations for the dependent variables at each node in the mesh. The dependent variables include the nodal masses of all five species, the unconstrained (stress-free) nodal volume expansion rate, the z and y components of nodal velocity, and the z and y nodal displacements. Equations for nodal masses and volume expansion rate are obtained from a control volume finite elementdiscretizationof the mass conservation equations (see Section 5.1). Equations for nodal velocities are obtained from a traditional Galerkin one-point-quadrature finite element discretizationof the momentum equation (see Section 5.2). Typically, nodes are sharedby several adjacent elements. Thus the elementcontributionsto a nodal variable equation must be superimposedin order to determinethe complete equation for each nodal variable. Adding these contributions is commonly referred to in the fite element literature (e.g. [34]) as the “assembly process.” After the assemblyprocess is compete, explicit time integrationfor the nodal masses,unconstrained nodal volume rate, and nodal velocities can then take place. The stress rate in each element is computed at the single Gauss-point horn the deformation rate tensor and the Gauss-point interpolated unconstrained volume expamion rate (see Section 5.3). The stress, z displacement, and g displacementare obtained horn a straightforwardexplicit time integration of the stress rate, z velocity and y velocity. Explicit time integration and time stepping is discussedin Section 5.4. Before proceeding to Sections5.1 md 5.2, it is usefulto introduce a few concepts. The four node quadrilateral elements used in the numerical implementation of the silylation model employ the following equation for spatial interpolation: 33 4 r? = ~ (54) QINZ 1=1 where 1 is the node number, IV is the bilinearshape function, and @ is the nodal dependent variable (mass, volume rate, or velocity). The global z – y coordinate system shown in Figure 9 is transformed into the isoparametric local ( - q coordinate system shown in F@ure 10. Mapping of points between the z —y and ~ —q coordinate systems is accomplished using the following equations: (55) (56) where the btiear shape function NI is and Z1 and 91 are the x and y coordinates of node I. The values ~1and ql are the coordinates of node I in the f - q pkme shownin Figure 10. Dependingon the node number, these values are either –1 or 1. 5.1 Mass !Ih.nsfer and Chemistry The species mass conservation Equations (10)-(14) and the unconstrained silylation volume expansion Equation (15) described in Section 4.1 have been discretized and solved on lD and 2D domains using two diflerent numerical methods. In the absence of stress, Equations (10)-(15) have been solved on a lD domain using an implicit ODE integrator, DVODE [35],which is a double preckion variable coefficient ODE solver for stiff or nonstiff problems. The solution (not shown) compares well with the solution of the model formulation described in Appendix A. The 2D formulation described here has been solved on a ID domain in the absence of stress; again the resultsagree well (cf. Section 6) with those obtained horn the solution of the ID formulation. Equations (10)-(15) have been discretized on a 2D domain using the control volume iinite element method (cf. [36], [37’J).This method has been applied to ablation problems with moving control volumes [38]. The control volume finite element metlmd (CVFEM) applies 1A conservation principlesto the finite element method. Figure 11 shows a four node quadrilateralelement that is divided into four regions labeled SCV1 through SCV~ these regions are associated with local nodes 1 through 4 of the element. Each of these nodes 34 corresponds to point P of the control volume shown in Figure 4 (note that if the node is on a boundsry, it will not be at the center of the control volume); e.g., if node 1 corresponds to point P then SCV1 is part of the control volume VP(t) associated with node P, and other parts of the control volume VP(t) lie in adjacent elements (not shown in Figure 11). The conservation Equations (10)-(15) are discretized on the control volume; however the matrix is assembled on an element by element basis as done in linite element method applications. As an example, the element matrix equation for the conservationof mass of unsilylated polymer resistmaterial C can be written as [Mc]{ac} = {Z)c} (58) where ?&l {6.} = ~ (59) 3 {} ?& 4 (60) and (&l V&q WC2 v&V2 {b.} = ticsv&?3 {} “ (61) WC4 V-4 Although the left hand side of Equation (58) is written in element matrix form, it is not sssembled into a global matrix because the solution method is explicit and the left hand side of Equation (10) is simply the time rate of change of mass of C. The right hand side of the elementmatrix equation, Equation (61), is assembledfrom the element level into a global matrix for the mass of C at the nodes by summing the contributions of all the elementsto each node. Siiar expressionsapply for material species D and 1? . For the mass conservation Equations (13) and (14) of species A and 1? and for the unconstrained volume Equation (15), fluxes across control volume boundaries must be evaluated. Referring again to the element shown in F@re 11, the surface areas of the control voh.rnesacross which fluxes are determined are labeled ssl through ss~ the midpoints of these lines are designatedas integrationpoints ipl through ip4; the fluxes are evaluated at the integration points. For example, for the subcontrol volume SC!Vl of Figure 11, the contribution of the mass difFusionflux, given in Equation (21), of species A to the total mass diffusionflux of speciesA (for the controlvolumecenteredat node 1) is 35 j. . ndA = ,,1j. - mid – j. . ndA I Scvl / JSS4 (62) where j. “ndA = jAzliPIAYl – JSsl jAvlipI& (63) and / .%4 h Equations (63) and (64), Azl j~ - ndA = jA=li@AyA – jAVlip4AxA = z= – ZmPl, A~l = ~C – !hr@, Lz4 = % and A~A = Y=– y~PA, where the coordinates of the points (x=,vC), (~rnpl,~rnpl), and (~mp4, (64) – %n~, gmp.4) are determined using Equations (55) and (56). The mass diffusion flux components of species A at the integration points are then given by (cf. Equation (21)): jA = ipl jA v ipl jA = –pDOetwpEiME1 ipl = –@e(WPE/ME) = –P~oe(wpE (66) ipl /ME) = @ jAv (65) (67) @ = –_p~oe(WE/ME) (68) @ w where the shape function derivativesare given by: [g]=[r[g] .-I+$ ;] 36 (69) (70) and ]JI is the determinantof the Jacobian of the transformation,given by Equation in the next section. Similarexpressionshold for the massdiffuion fluxcomponentsof B at integrationpoints ipl and ip4 as well as for the mass -Ion flux components of A and B at integration points ip2 and ip3. The area integrals and source terms on the right hand sides of Equations (13), (14), and (15) are thus evaluated, yielding a right hand side vector for each element. The element right hand side vector is then summed over all nodes in the problem to yield the global right hand side vector, which is the discretizationof the right hand side of Equations (13), (14), and (15) for each control volume in the problem. The discretized equationa for the conservation of mass of species A, B, C, D, zmd E at each control volume are then integrated in time as discussed in section 5.4. The discretized equation for the unconstrained volume is not integrated in tim~ rather it is the time rate of change of the unconstrainedvolume (the discretized right hand side of Equation) 5.2 that is required by the materialmodel and discussedfurther in section 5.3. Momentum Chain rule di.flerentiationof the inertia term in Equation (27) yields ]V+p+=v”a+l?. (71) As discussed in Section 4.2, the material motion is quasi-static (negligible inertia). However, to facilitate the time integration for velocity, v, only the first term in Equation (71) is neglected, resulting in P+= V.a+F. (72) The method of weighted residuals (see e.g. [34]) is applied to Equation (72) in order to obtain expressions for the contribution of each element’s the nodal velocity components. The process begins by casting the momentum equation into the so-called “weak-form” by multiplyingits residualby a set of weightingfunctions, W, and integrating the result over the volume of an element, i.e., /v (p+- F- V.c)WdV=O. (73) The method of weighted residuals seeks to approximate the dependent variables (nodal velocities) using simple functional forms like that expressedin Equation (54). The functions are selected in such a way as to minimize the residual in Equation (73) thus approximating the solution over the element domain. For the planar problem solved here, a special form of the method of weighted residuals rekrred h as Galerkin’s method is used. Jn Galerkin’s method, the weighting functions are chosen to be the same as 37 the shape or interpolation functions in Equation (54), i.e. W = N. Furthermore, for an assumed element thicknessof unity, cW = dA, where A is the element area. Thus for the planar problem: (p+- F- V.a)NdA=O (74) /A Integration of the first two terms in the above equation is straighforward,but the last term requiresspecial treatment because of the spatial gradient operator V. The last term is integrated by parts using Green’s theorem. This resultsin two integralsfor the third term: m area integral and a line integral, where the liie representsthe perimeter of the element, i.e.: PNdA-LFNu+l””vN’q”” fids’o (75) The line integralsprovide a means of applying the appropriate global boundary conditions. They vanishfor all elementswhich do not have a physical boundary (i.e., all intenor elements). The vector n representsthe outward pointing normal to the line S (element side on a physical boundary). Hence in Equation (75), the last term representsthe iniluence of boundary traction and pressureon the boundary elements. As discussed in Section 4.2, this term is zero for the elementsalong the free surface for the silylation problem. Integration of Equation (75) proceeds in the & - q coordinate system (Figure 10 aad Equations (55) through (57)) by noting that dA = dxd~ = IJldfdq. (76) where IJI is the determinant of the Jacobian transformationmatrix, i.e., (77) It cau be shownthat the determinant of the Jacobian is directly related to the elementarea, IJI= ; (78) where the element area is expressed in terms of the elementnodal coordinates as (79) Substituting ~ - q transformations into Equation (75) and neglecting the boundary integralsyields, 38 +1 +1 +1 +1 +1 +1 pvN IJIo!@q – a . VN IJIqdq FN IJId~dq + /1-1 –1 // -1 –1 –1 11 -1 = O. (80) Integration of the terms in the above equations is performed numerically at a single Gauss point for which &= q =0. The ‘one point quadrature” assumptionis a popular simplification (see e.g. [16, 39,40, 41]) that permits p, F and a to be taken as constant over the elementand moved outside of the integral. The element stress tensor, a, is the time integrated value of the elementstress rate tensor which is determined from the deformation rate tensor and the unconstrainedvolumetric strain rate. The element deformation rate tensor is, in turn, computed from the nodal velocities. The remainingfunctions Wide the integrals such as ~1 are easily integrated since, e.g., Nz(g,q) = NI(O,O) = ;. (81) where I=1234 >>). After some algebra, Equation (80) takes the form, ~[~{+} = ~{F} + :[AI{o}. (82) The mass matrix [M] is given by [M] = 10101010 01010101 10101010 01010101 10101010 01010101 (83) 10101010 01010101 10101010 01010101 The “lumped” matrix approximationis made for the mass matrix, i.e., each row of the mass matrix is summed and the result placed at the diagonal. Hence, [Ml = 4[q where [~ is the identity matrix. ‘ The column vectors {+}, {F} and {a} and the matrix [A] iu Equation (82) are given by 39 (84) U1 u VI U2 {+} = (85) g V3 ii4 V4 (86) rs== {CT}= (87) ~ {} O*: ~2–Y4 0 x4–x2 O- X4 – X2 O ~2–’#4 o ZI— Z3 O o ?J3-Y1 [A] = 0 y4–!/2 0 m-w 0 xl –X3 ‘y3-gl o O Z2–Z4 O X2 ‘x4 o x3 — q O 34–Y2 Z3– (88) Z1 l/1-?J3 o 0. Equations (82) through (88) represent a set of eight ordinary diilerential equations that determine the contribution of each element to the nodal accelerations UI and til where 1 = 1,2,3,4. Once the total contribution is determined horn the fm.iteelement assembly process.,Ex and 51 are integrated with respect to time to yield the nodal velocities tiI and tiI. The nodal velociti~ are then integrated with respect to time to yield the nodal displacements Uz and vI. For the silylation problem the body forces Z?zamdFV due to gravity are neglected. SpatiaJintegration of elements at a single Gauss point (i.e. the one point quadrature used here), though computationally efficient, can lead to the calculation of spurious energy modes. If these modes are not supressed, a severe and usually fatal mesh distortion lmown as “hourglassiig” may occur. The term hourglassingwas coined because of the shape taken by adjacent tiected elemmts in structuralcalculations. The problem arises when a single gauss point is used to calculate the deviatoric part of the stess. When this is done, it is possible for certain element deformations to produce no stress. Such deformations are commonly refered to as ‘zero energy modes.” 40 —. .—— ..—.. —..—— A number of procedures have been developed to eliminate hourglassing (see e.g. [42, 43, 44]). In the present work, the method of Flanagan and Belytschko [43, 44] was used. Details of their procedure and its implementationinto the silylation model will not be presentedhere. 5.3 Polymer Material Response AS is typical for constitutive laws involving inelastic deformation with a evolving microstructure, the equations for the material point stress and state variables form a coupled system and must be solved simultaneously.To do this, both Equations (42) and (47) are discretizedin time and coupled using Equation (43) for the plastic strain rate. 5.3.1 Coupled Newton-Raphson Algorithm A hid stress state is computed by assumingthat the strain step is entirely elsstic. Using (37) and (38) Itrial 0 =c’n+C” -Dp : Q’ -{cz~+&}qAt (89) and U’n+l= u‘t’id – ~pc” : NAt (90) Use of the associativeflowrule allows(90) to be reducedto an equation for the stress magnitude alone in terms of the as yet unknown plastic strain increment (91) Jn order to satisfy the consistency condition for yielding at the end of the time increment, horn (43) ,p=,oq[-Ayl{l-(-)5}] (92) Finally, the strain softening state variable evolution depends upon the pk@c strain increment (93) 41 Now Equations (91) and (93), coupled through the plastic strain rate magnitude, are formed into residuals usinga Newton-Raphson iteration scheme which, after whkh must ultimately vanish. This is accoxnplished some manipulation, results in the matrix equation (94) where (95) (96) (97) (98) are the components of the Jacobian matrix and ‘2’=’Oq[-A~l{l-(-)’}1 (99) is the plastic strainrate iterate. Equation (94) is solvediterativelyfor the incrementsin stress,flow resistance aud plastic strain. 5.3.2 Radial Return Algorithm While the Newton-llaphson algorithm is generallyvery fast, it can e.xhibkd.i.f%culties when the local Jacobian vanishes. For strain softening polymers, this can happen during the @lc-pWlc transition where the local material stifbwss tends to zero. At the expense of somewhat less accurate simpl@ng assumptionswhich 42 vanishwith decreasing time step size, a simplerradial return algorithm does not suikr this drawback. While a general radial return algorithm would necessitate iteration owing to the nonlinear constitutive law, we examine a simpler case for illustration here. For glassy polymers below the glass transition temperature, ~ ~ 1 s. ~ t~ ~it, a s~gle iterate r~ retm is possible. TO fo~tiate this algorithm, however, a quasi-yieldfunction must be deiined by inverting the flow rule. Writing (43) in terms of the stressmagnitude required to produce strain rate 5P for ~ = 1, ‘-{s+ w$))=o (loo) Now (100) takes on the form of a yield function Y=T–R (101) where the yield resistance, R, has two contributions: a part, s due to flow resistanceinherentin the material structure and a rate dependent contribution (102) Taking an elastic strain step4 (103) ‘trinl _ – Sn – h~lDIAt (104) The yield condition (101) is then checked. If Y <0, the step is elastic and the solution triat Tn+l = 7- Sn+l= Sn (105) (106) 3For the representativepolymer PMMA, for irrstzmce,the ratio ~ = ~ . 41t has been shown [32] that, for large in-lc strain problems of this type, it is a good approximation to replace 7P with [Dl in the yield condition and the state variable recovery terms. 43 is accepted. If Y >0, then the consistency condition Y = Omust be enforced in order for the stresssolution to remain on the inverted yield surface. This results in an expression for the appropriate plastic strain increment +11 &P = ~PAt _ &ial = 2p+h -b (107) which is then used to update the stress and the flow strength Tn+l trial = 7- – 2/..LAyp Sn+l = Sn+ hA~p 44 (108) (109) 5.4 Time Integration Solving for the dependent variablesin the silylation model is accomplished by integrating first order diEerential equations in time. The general form for each of these equations is g= RHs. “ (110) where @ representsGauss point stress or nodal velocity, mass, or unconstrainedvolumetric expansion rate. The shorthand notation “HIS” is used to indicate the right hand side of the equation which in general is a function of other dependent variables. A simple “upwind” difference in time is used to approximate the time derivative, i.e., @n+l gya _ = RHS. (111) At where n + 1 refers to the value at the new time step and n the previous time step. Whether the integration is explicit or implicit depends, of course, on how RHS is evaluated. If it is evaluated at n +1, then the integration is implicit and a nonlinear algebra problem must be solved. If it is evaluated at n, as in the present model, the integration is explicit, resultingin the following simple expression for @n+l, Q“+l = @n + (At)RH7 (112) Explicit integrations are conditionally stable which requires careful selection of the time step At. Courant stabtity limits for the mass transport (-Ion) The and mechanical calculations are functions of element geometry (shape and size) and the properties of the material. For quadrilateralelementsit cm be shown that (113) @2 At~ = T where& (A+ (114) 2/L)B and At~ are the stable diflusionand mechanicaltime steps,respectively.In the aboveequations‘D is the diffusion coefEcient,A and p are Lam@=Ic constants, A is the elementarea (given by Equation (79)) and B is an element geometric parameter given by, ~ = (?42– !/4)2 -1-(?/1– !/3)2+ 2 45 (21 – Z3)2 + (Z 2 – 3 4 )2 (115) -. Generally Atd is orders of magnitude larger than Atm. Hence, the time step is controlled by the mechanical part of the problem. In quasi-static mechanical calculations such as this one, an =bitrary scaling of the mechanicaldensity can be used to bring the stable mechanicaltime step more in line with the stable dfilon time step. Hence for the purpose of computing material momentum, the density can be increased by some arbitrary scaling factor, ~d, where jd >>1. The effect of this density scaling is to add artficial inertia to the momentumequation. Care must be taken to insurethat density scaling does not compromisethe quasi-static behavior of the momentum equation. When this occurs, the solution exhibits obvious nonphysical behavior. In order to avoid excessivelylong computational times, some measureof mechanicaldensity scaling was used for all silylation calculations presented in Section 6. 5.5 Mesh Generation and Post Processing A Fortran77 computer program was written to solve the 2D silylatiotiequationspresentedhere. The program reads a keyworddriven problem descriptiontie usingthe parsing languagedeveloped by Perano and Kaliakin (see e.g. [45, 46]) and generates a two-dimensional computational mesh using the Sandia developed mesh generator AP [47]. During the solution phase a seriesof time-state files are written, typically one time-state file for each second of simulatedtime. Post processing of these time state files is also performed using AP. This relatively mature pre and post processing interface makes it possible to set up, solve and post process a wide variety of 2D silyation problems quickly and accurately. 46 6 Results Calculations made using the previously described models are presented in this section. The first two subsections deal with model verification for simple geometries including transientbehavior at a single point and transient behavoir in one spatial dimension. The remaining subsection discusses the fully coupled transient behavior of the silylation model (2D Model) in two spatial dimensions. 6.1 Constitutive Model Material Point Simulations The polymer constitutive model consists of the coupled system of Equations (42), (43), and (47). They are parametrizedby a set of materialconstants {p, A,TO,sO,sSS,A, h, p, q}. These parametersare generallyeither derived from other physical quantitieswhose values~e known for a classof materials(e.g. see Equation (44) for A) or fit to experimentaldata. A complete set of parametershave been fit for the amorphous polymer, polymethylmethacrylate (PMMA) elsewhere [30]. These are used here to illustrate features of the model behavior for a material point under three distinct homogeneous deformation paths. We exercise the model for constant strain rate uniaxialextension to show that it recreatesthe intermediatestrain tensile data. We then use the model to simulate one-dimensional free surface swelling cha.rzwteristicof the behavior of the silylating layer far from the substrate edge. Lastly, we artificiallyrestrict the model to elastic deformations only to illustratethe effect of non-volume conserving flow on the degree of free surface swelling. For all the discussionswhich follow, the material properties correspond to a constant temperature of @ = 90 ‘C e Q9 and are given by /.L=lxlo10 A=3.33X109 30 = So = dyne ---&- ~ (116) (117) x 1011 S–l (118) dyne 8.8 X 10s ~ (119) 1.0 (120) 47 A = 1.7X 10–5 h=9x10g P=!l=~ 6.1.1 K–cm2 dyne (121) ~ (122) (123) Uniaxial Extension The viscous inelastic softening response prior to polymer locking behavior can be paraneterized entirely by scalar state variables in the model. These can be estimated by fitting the model to uniaxial deformation paths. Here, we adopt the set of material parameters in [30] listed above. These have been chosen, for illustrativepurposes, to reproduce the uniaxial tensile data for PMMA [30]. The response of PMMA pulled in uniaxialextension at a rate of unity is shown in F@re 12 along with the responsereported by Boyce et al. [30]. Only the softening portion of the response is included in the model presentedhere. Fbrther, the ratedependent yield and post-yield softening response are functions of the applied deformation rate magnitude IDI, the hardening modulus h, and the steady state inelastic flow strength S8,. The mechanical response depends on these model parametersas illustratedin Figures 13,14, and 15 respectively. A similarlyisolated integration of the model equations is compared in FQure 16 with its implementationin the twmiirnensional finite element code. 6.1.2 Free Surface Swelling Because we interihce the material model with a displacement-basedfinite element code, the constitutive subroutineis generally driven by a strainhistory discretizedin time. This is straightforwardwhen kinematic or Dirichlet boundary conditions are imposed on the finite element mesh. When Neumann, flux or stress boundary conditions are imposed, iteration is generally necessary in order to converge on the appropriate strain step necessary to satL@ this boundary condition. This is the case for the onedimensional silylation experiment described in Appendix A where the free surface of the silylating layer is traction-free. In order to test the constitutive algorithm for this case, an iterative loop was added which imposed a stress free boundary condition in one-dimension. To do this, an elastic strain step is presumed for the out-of-plane swellingresponse (124) This strain step is applied, resulting in a non-zero normal component of the stress at the layer surface. The strain step is iteratively adjusted until the normal stress component vanishes. At each increment, the corresponding compressivein-plane layer stressesaccompanying the induced silylation swelling are given by AOII = AU33 = –2pAe22 (125) The deviatoricresultant in-planecompressivestress associatedwith a constant rate of swelling,~, is shown in Figure 17. Note that the response is similarqualitativelyto the cinematically driven response describedin the previous Section. After yielding, the polymer inelasticflow strengthevolves in accordance with Equation (47). This state variable evolution is depicted in Figure 18 for one-dimensionalfree surface swelling. 6.1.3 Idealized Elastic Swelling As a final case study, we restrict the mechanical constitutive model to purely elastic behavior. For this example, the free surface expansion can be computed in closed form as a means of verifying the numerical implementation. Such a case will only be physically plausible for very small deformations. To simulatethis behavior, however, the yield strength, .s0, can be artificially increased to ensure that no plasticity occurs. For this limiting case, the incrementalrelations (124) and (125) can be replaced by total quantities and the free surhce strain will be given by (126) and the accompanying lateral compressive stressesin the silylation layer are (127) both in agreementwith the one-dimensionalmodel presentedin Appendix A. The dependence of free surface swellingon Poisson’s ratio is shown in F@ure 19. For swellingsto large strain, the fully elastic-plasticfree surface motion is closely approximated by the elastic limit as v + ~. The isolated model integration is shown in this limit in [email protected] 20 along with the twdirnensional finite element model implementation. 6.2 Comparison of ID and 2D Model Results A simple one-dimensional transient problem was formulated for the purpose of comparing the 1-D modd discussedin Appendix A and the more general2D silylationmodel. It is assumedthat a layer of unsilylated 49 and uncrosslinkedresist 1 x 10–5 cm thick is applied to a perfectly rigid surface. The resistand the surface it is atthed to are tite in the plaae orthogonal to the thicknessdirection. At t=O a fied concentration of the silylation agent is imposed on the free surface of the resist. As time proceeds the agent diflkses into the resistfrom the free su.rfcweand reacts with the resist, causingsilicon deposition and swellingof the entire resist layer. Material parametersand boundary conditions are selected such that an unstressedlayer would double in thickness when the layer becomes fully silylated. Figure 21 illustrates the mesh and application of boundary conditions for the transient 1-D problem. A single row of 100 elements was used in the 2D simulation. (Correspondingly, 100 grid points were used in the model described in Appendix A). The lower two nodes of the bottom element (at z = L(t)) are fixed to the rigid surface and no mass tramsferis permitted across the elementface. Symmetry boundary conditions are imposed on the right and left boundaries of the element row to simulatethe “i.dnite layer”. Constant values for the concentrations of species A and E are imposed on thefree surface. The coordinate z shown in the figure is the nondimensionaldistance (i.e., actual distance divided by current length L(t)) from the free end. It should be noted that the model of Appendix A does not account for the presence of a silylation reaction product speciesB. Furthermore,the constitutivemodel used in Appendix A relatingstressto strainis a simple elastic model and therefore not realisticfor real silylation problems undergoing large strains. Nevertheless, the model is useful for ver@ing coupled momentum/dif7usionbehavior in ID. Dimensional parametersused in the 2D model are related to the dimensionlessparametersof the lD model according to the following: k.p. T (128) ‘Et p_ PA (129) PA= (130) (131) pz MC B=—— ~C M. 50 (132) (133) (134) r~ = k8pA ~& = 1 .() (135) A (136) where, Pc — = 0.0113207mol/cm3 co= MC (137) Note that if the silylation reaction rate does not depend on stress, then b = O in Equation (180) and kl = k. = ks. Also note that the d.ifl&ion coefficient, D, in Appendix A is equal to DO in Equations (21) and (22). In Equations (128) through (136), ~ is dimensionlesstime, and P, Q, ~, and ~ are dimensionlessdependent variables. The remainingdimensionlessparametersdepend on materialproperties and boundary conditions and remain constant throughout the calculation. The dimensionalproperties used in the 2D calculation were determinedfrom the constantsspecifkd in Equations (128) through (136) and can be summarizedas follows: pA = 39.62 g/cm3 (138) Pc = 1.2 g/cm3 (139) p. = 2.6717 g/cm3 (140) p.= (141) 1.3359g/cm3 51 MA = 175 g/mol (142) M= = 106 g/mol (143) M. = 236 g/mol (144) M.= (145) 236 g/mol p.= = 19810 g/cm3 (146) k,= 175 cm3/(mol –s) (147) DO= 2 X 10-12 cm2/s (148) LO= 1 x 10–5 cm (149) ke = 1.981 S–l (150) Comparisons of predictions horn the ID and 2D models are shown in Figures 22 ad 23. Figure 22 shows the growth of the silylation layer as a function of time for two diflerentPoisson ratios, v = 0.3 and v = 0.5. It can be shown that the final layer thicknessescorresponding to v = .3 and LJ= .5 are 1.62 x 10-5 cm and 2.00 x 10–5 cm respectively. The latter value follows from the observation that all elastic deformation for the case where v = .5 must be volume conseming. In each case, the 2D model correctly predicts the ihal layer thickness. The ID model of Appendix A predicts slightly lower values of layer thic.lmessthan the 2D model throughout the transient. These di.ikrencestend to disappear when more grid points are used in the lD model. The dashed blue line shown in Figure 22 illustratesthe effect of including plasticity in the material response. The 2D model prediction of the layer thickness for v = .3 shows that when plastkity is included, the growth of the silylation layer is nearly vohrme conserving. This is the case because only 52 a small portion of the total deformation is elastic. The volume-preservingplastic deformation comprises a substantiallylarger percentage of the total deformation. Figure 23 compares ID and 2D model predictions for mass concentrations and stress distributionsat three diflerent times during “elastic” silylation with v = .5. Figure 23 (a) shows the mass concentration of the silylation agent at 2.5,10, and 60 seconds. Results ue plotted using the nondimensionallength z(t). Siar comparisons for the mass concentration of unsilylated resist and silylated resist are shown in (b) and (c). In A case the comparisons are excellent. The lateral stress distributions (a== and C=Z) are compared in Figure 23 (d). Each model predictsa finalcompressivelateralstressof 1.05x 1010dynes/cm2. Transientstress distributionsalso compare favorably, although the 2D model predicts a small compressivestressthrough the layer thicknessduring the propagation of the silylation front. 6.3 Transient 2D Silylation - Nominal Calculation In this subsection the results of a 2D transient silylation simulation will be presented. The objective is to simulatethe materialbehavior observed during a process similarto that shown in Figure 2. Figure 24 shows the computational mesh for the nominal 2D silylation problem. A resist4 x 10–5 cm in thicknessis bonded to a rigid surface (fixed surface). Part of the free surface has been crosslinkedto a depth of 6 x 10-6 cm. The planes of symmetry shown in Figure 24 are used to simulate a series of silylation lines parallel to the z axis. The distance between any two adjacent silylation centerlinesis 4.8 x 10-s cm. In this simulation, crosslinkedbehavior is modeled by reducing the di&usivity,’DO, of the silylationagent in crosslinkedelements by three orders of magnitude relative to the uncrosslinkedelements. Except for the dif?usioncoefficient, all other crosslinkedand uncrosslinkedmaterial properties are assumed to be identical. These properties are summarized as follows: E = 2.25x 1010dyne/cm2 (151) v = .125 (152) 10s dyne/cm2 (153) h = 9.0x 109 dyne/cm2 (154) A = 1.7x 10-5 K – cm2/dyne (155) sSS= 7.7x 108 dyne/cm2 (156) c1= 0.0 (157) to = 1.0x 1011s–l (158) M. = 175. g/mol (159) so = 8.8x 53 g/mol (160) 34= = 106. g/mol (161) MD= 236. g/mol (162) M. = 236. g/mol (163) ~.4= 39.62 g/cm3 (164) pB = 100.0 g/cm3 (165) p= = 1.2 g/cm3 (166) p~ = 1.2 g/cm3 (167) p. = 1.2 gfcms (168) k, = 2.0 S-l (169) k. = 175 cm3/(mol –s) (170) w = Ocm3/mol (171) CrosslinkedDO= 2.0 x 10–15 crn2/s (172) UncrosslinkedDO= 2.0 x 10–12 cm2/ s (173) A4B = 45. Initiallythe resist is assumedto be unsilylated; i.e., it is composed entirelyof species C. At time zero the fiwe surface is exposed to the silylation agent, species A, and mass transport, chemistry,and material swelling proceed for a simulated time of 60 seconds. Boundary conditions for mass tramfer may be summarizedas follows: .4t the free StiaCe: CA= .05, ~E= .95 At the right s~etry boundary No mass transfer At the left sprnetry boundary No mass transfer At the bottom fked boundary Nomass transfer Boundary conditions for momentum may be summarizedas follows: .4t the free surface Stressfree, see Equation (28). At the right symmetry boundary: d = O At the left symmet~ boundary: d = O 54 At the bottom fixed boundary: u = 6 = O F@we 25 shows cross-sections of silylation agent (species A) concentration and material swelling for time states corresponding to O, 1, 3, 5, 10, 30 md 60 seconds. To aid in visualization,resultshave been mirrored across the left s-etry boundary which corresponds to the y axis and the centerlineof a typical silylation cross-section. At t=O seconds, the resist is unsilylated except for the imposed concentration of A at the free surface. As early as t=l second, evidence of species A -Ion, silylation reaction, and the onset of material swelling is clearly vk+ble in the uncrosslinkedregion. As time proceeds the silylation agent continues to diifuse preferentiallyinto the uncrosslinkedresist causing a localized material swelling similar to that observed in Figure 2. In this particular simulation, the effect of crosslinkingwas modeled using a reduced diffusion coefficient for the crosslinkedmaterial. As a result, masstransport, chemistryand material swelling in the crosslinked region were localized to a very thin layer along the free surface. This layer is most visible at t=60 seconds and resemblesthe swa,lled “scum” layer observed over crosslinkedregions after silylation experiments. Figure 26 shows concentration profiles of the silylation reaction product (species B), the unsilylated solid (species C), the unexpanded silylatedsolid (speciesD) and the expmded silylatedsolid (speciesE) at the end of the simulation (t=60 seconds). As expected the highest concentration of the silylation reaction product (species B) occurs above the reaction front. The contour plot of unsilylatedsolid (species C) clearly shows the depletion of species C above the reaction front. The blue core representinglow concentrations of C shows that the reaction front is beginning to move under the crosslinkedportion of the resist. Because of the relatively short relaxation time for the conversion of species D to E (k, = 2.0 s–l ), significantamounts of D are confined to a relatively narrow band (shown in green in the D concmtration plot) near the reaction front. The portions of the resist which have completed the conversion to the expanded silylated solid are clearly visible as a yellow band in the concentration plot for species E. Stresses developed during the silylation process ultimately determine the amount of swelling and the final shape of the material. Figure 27 shows cross-sectionalcontours of stressat 60 seconds. Included are contours for the through-thicknessstress,avv, the lateralstress,cr.=, and the shearstress,cr=g.Note that the individual in-plane deviatoric stress mmponent, cr~z,reaches values beyond the rate-dependent yield strength. This indicates that substantial inelastic flow can occur at swelling rates typically seen in the silylation process. We point out for clarMcation that this is true if the polymer plastic properties are reasonably represented by those for PMMA. When there is substantialinelasticflow, volume conservation dictates that the amount of swelling will be only marginally affected by the stress response. The stressesand ldnematic constraints, i.e. the geometry of the crosslinkedregion, however, can play a substantialrole in the overall shape change of the material. The relationshipbetween crosslinkedand uncmsdinked regionsat the conclusionof the simulationis shownin 5s FQ_ure28. Elementsof the crosslinkedpolymer are distinguishedsolely by their reduced difFusioncoefficient relativeto the uncrosslinkedmaterial. Elementsin the crossli,nkedregion are shown in blue and elementsin the uncrosslinkedregion are shown in red. The original shapes of these regions prior to silylation are shown in Figure 24. As expected, nearly all material swellinghas been confined to the uncrosslinkedregion. Some swellingis evident in the crosslinkedareaalong the free surface (the “scum layer”) and in the intenor interface between the crosslinkedand uncrosslinkedregions. This was to be expected since, in this modeling scenario, the crosslinkedregion was assigned a finite (but low) Hlon coefficient and the silylation reaction was permitted to proceed unimpeded when the si,lylationagent was present. Here, the crosslinkedpolymer has been assigned the same constitutive parametersas the uncrosslinkedpolymer. Figure 28 shows that in this case, silylationextends to the materialbelow the crosslinkedregion. The resultingswellingcausessubstantial upwardmovement of the crosslinkedmaterial. The simulationsuggeststhat such movement is likely to occur whenevercrosslinkingis suilicientlyshallow and silylationtimes are long. It may be expected, however, that these responsefeatures would be sensitiveto changesin the inelasticmechanical properties of the crosslinked region. For future reference, it may prove worthwhileto implement the locking Langevin spring element in the large strain model (this may representcrosdi.nkedbehavior as described in Section 4.3.5) or perform a parameterstudy over a range of stii7ereffective elastic modti for the crosslinkedpolymer. 56 7 Concluding Remarks This report documents the developmentof a new two-dimensionaltransient finite element model of the positive tone silylation process. The 2D model focuses on the part of the process in which crosslinkedand uncrosslinkedresist is exposed to a diflusiig silylation agent. The model accounts for the combined effects of mass tram.port (silylation agent md reaction product), the chemical reaction leading to the uptake of silicon and material swelling, the generation of stresses, and the resulting material motion. Both Fickian and Case II Hlon models are incorporated. The model provides for the appropriate mass transport and momentum boundary conditions and couplesthe behavior of uncrosslinkedand crosslinkedmaterialsas well as the compliance of any underlying device topology. The 2D model behavior has been veriiiedusingindependent materialpoint simulationsand a one-dimensional transient finite difference model. Finite element mesh generation, problem setup, and post processing of computed results is sufficiently mature to permit the 2D model to be used in the investigation of a broad parameter space which includes materialproperties and geometry issues (e.g., alternate crosslinking distributions, resist thickness,etc.). The authors recommend that future work in the area of silylation modeling be directed toward exercising the present model over a range of parameterspace to determinewhich properties and phenomena have the largest affects on the silylation process. In addition, studies should be undertaken to determine relevant material properties. In order for the model to become a predictive design tool, it is recommended that it be verified against a set of well controlled silylationexperiments. 57/58 A APPENDIX - A One-Dimensional Model A logical first step in simulating a process such as silylation is to develop a one-dimensional model. This can be used to obtain a basic understanding of at least some aspects of the process with a minimum of computational effort, and it can also serve as a validationtool for multidimensionalcodes that will eventually be used to describe the process in detail. The model to be presented here is similar to that developed by Zuniga and Neureuther [19], but there are some difFerenc~ in implementation. In partitiar, the bulk (convective) flow arisiig from the volume change of the solid is now accounted for explicitly in the species conservation equations, and simple analytical relationsfor the amount of swelling and the induced stresses are derived. However, since the stress-strainrelationsare based on the linearizedequations of elasticity,the model is strictly valid only for small deformations. If it is assumed that the gas phase is well mixed aud nonreacting, then the modeling effort can be confined to the solid, According to Zuniga and Neureuther,there are four chemicalspecies of interest: the (dissolved) silylating agent P, the unreactedhydroxyl group Q, the silylated but unexpanded group U, and the expanded group E. In what follows, the latter three specieswill be assumedto include the associatedpart of the polper chain, so that all of the solid material is accounted for. Clearly, only species P is bee to difhse, and it can react with Q according to P+ Q~U+... Here kl is a second-order rate constant, provided that mass-action kinetics is appropriate. While there may be mobde reaction products in addition to the fixed group U, they will not be important as long ss their volumes are negligible md the reaction is irreversible. In any case, the newly formed U is eventually converted to E with a characteristicrelaxation time tr: The primary task is to determine the tim~dependent concentration profiles for the four species and, of particular interest, the total amount of swellingas a function of time. The differentialequationsgoverningthe process can be obtained by applying the generalmaterialbalance [48] to each species. In molar units and in one dimension, (174) where t is time, z is the spatial coordinate (measuredperpendicularto the gas-solid interface), v is the massaverage velocity, Q is the molar concentration of species i, Ji is its diilusion flux, and ~ is its volumetric production rate by chemical reactions. Letting CP- P, and so forth, one has 61P ~ ~ap at —= ax _pav—— aJp — – klPQ 8X ax 59 (175) (176) (177) au au T#vz= au u -“% +k~pQ–z (178) Following Zuniga and Neureuther,we approximate the difksion flux by Jp = –D ~wE ~ (179) where w is a parameter that allows the dif7usivity of species P to vary with the extent of swelling, D is simply the -lvity in the initial state. In addition, at fixed temperature the rate constant kl is assumed to vary with the local stress a according to kl = b eb” (180) Equation (175) then becomes (181) and there are similar modifications to Equations (177) and (178). Since there are now four differentialequations but six unknowns (P, E, Q, U, v, and c), two additional equations are needed. These can be obtained by relating the stresslevel and the overall molar concentration of the f%iedspecies to the extent of the swellingreaction. The first step is to introduce a parametera, which is the ratio of the molar volumes of species E and U in a completely unconstrained system; it is assumed that the molar volumes of Q and U are identical Then, in the absence of stresses,the fractional increasein volume at any point (neglecting the effect of the dissolved gas P) is E(a – 1) ‘= ~ Q+U+E (182) The quantity@plays the same roleas a thermal expansionin the equationsofstatics. Forthe one-dimensional system being considered here, there is no stressin the z-direction, while there are no strainsin y and z. The equations of equilibrium [49] then show that the axial strain is C==@ l+V 3(1 – v) (183) and the transversecompressive stressis (184) 60 where Y is Young’s modulus and v is Poisson’s ratio. FkomEquations (182) and (183) it can be shown that the molar concentrations of the fixed solid species are related by Q+u+E(Y+E(a-l)y ::;)=@ (185) where @ is the initial molar concentration of Q. For the special case v = ~ this reduces to (186) Q+tJ+EcY=q which is identical to the result that would be obtained for a completely unconstrained system. Now, the initial conditions to be used with the differentialequations are Q=%, P= E= U=O, and v=O at t=O (187) Ltiwise, if the gas-solid interface is defied to be the position z = O,then two boundary conditions are P=P. and v=O at (188) z=O Here P. is assumed to be given by the volubility of the gaseous silylatingagent in the solid, i.e., equilibrium prevailsacross the interface. Analysis of the fit-order partial dii?erentialequations (176)–(178) shows that all of the associated characteristicbase curves originate on the line t = O in the t= plain+ therefore, the initial conditionson E, Q, and U are sufllicientto determinethese functionswithout any boundaryconditions being speciiied. However,a second boundary conditionis needed for P. If the thickness of the silylation layer is denoted by L, and if the solid substrate at z > L is impermeableto the difbing species P, then an appropriate condition is aP ~=0 at x=L (189) Finally, because the swellingrea&ion causes the silylation layer to thicken over time, L varies and must be determined as a function oft by means of the kinematic condition dL — = v(L, t) dt (190) together with the initial condition L=LO at t=O (191) The set of governingequationsand auxiliaryconditions is now complete. However,in order to assistin dealing with the moving boundary at z = L(t), it is convenient to introduce the new spatial variable z = z/L(t), so that the problem is now defined over the interval O < z <1. With this trtmsformation,Equation (181) becomes ~wE _8P =E(v-z~)=-~E+= and similarresultsare obtained from Equations (176)–(178). 61 () 82 – h eb”PQ (192) It can be seen that the mathematicalmodel involves eleven diilerentparameters,namely D, W, kO,b, t,, a, co, v, Y, p., and Lo. In order to reduce this number somewhat and to make the solutions more general, one can write the system in dimensionlessform using the following variables: The new equations need not be given here, but sufiice it to say that they involve seven dimensionless parameters (as could have been anticipated from the Buckingham Pl Theorem [50]): (194) All of the sample solutions to follow will be characterizedin terms of these dimensionlessgroups. The most straightforwardway to solve the system is via the method of lines. The first step is to replace all spatial derivatives with ii.nite differenc~, thii leads to a mixture of algebraic equations and timedependent ordinary di.flerentialequations. Upwind diflerencing is used in order to assure stability. The differential/algebraic system is then integrated in time using the packaged solver DASSL [51], which is ideally suited to problems of this kind. .4 sampleset of solutions will now be presentedfor an arbitrarilychosen set of pammeters. The quantities~ and ~ are both taken as unity, so that the diffusivity of the silylating agent increases with time while the rate constant for the silylation reaction decr~es. The values of the swellingfactor a and Poisson’s ratio v are 2 and 0.5, respectively; these should be fairly realisticfor the polymeric materialsused in this process. The concentration ratio ~ and the parameter T,, which is a ratio of characteristictimes for the swelling and silylation reactions, are set equal to unity. Finally, the parameter 6, which is a ratio of characteristic times for silylation and diflusion, is assigned a value of 0.01. The relativelyslow -Ion should lead to the appearance of a fairly well-defied reaction front that moves into the material as time progresses. The resultsof the computations are shown in Figures29-36. In all but the last two figures, spatial profilesare plotted for fixed moments in time separatedby AT = 0.5. Figure 29 shows that the concentration of P in the solid builds up only slowly,because the silylatingagent cannot &fFusevery far before being consumed by the relatively fast chemical reaction. This situation is alleviated only when the fixed reactant Q is consumed. Concentration profiles for the latter are shown in Figure 30, where the aforementioned reaction front is evident (although not particularly sharp). In Figure 31 one seestypical behavior for a reaction intermediate at any position, the concentration of U initially rises as a result of the silylation reaction, but it eventually erodes as U is converted to E. Since the concentration of U is generally small even at its peak, it is not surprisingthat the profles for E in F@re 32 are complementary to those for Q in Figure 30. However, it should be noted that the maximum concentration of E is 0.5 rather than 1.(), ss a result of the swelling factor a =2. The dimensionlessvelocity profilesare shown in Figure33. To reiterate,these nonzero velocities ariseentirely 62 as a result of the swellingreaction. Since the boundary at z = Ois fixed, the maximum velocity at any time must always occur at the other boundary z = 1 (z = L(t)), even if the reaction front has not yet reached that point. However,the value of the maximum velocity is a nonmonotonic function of tim.q this is obvious from the fact that it must start at zero and then eventuallyreturn to zero when the process is complete. F’or the case at hand, the velocity at z = 1 appeazs to reach a maximum at about I- = 2.5. The stress profiles in Figure 34 bear a strong resemblanceto the concentration profiles for species E in Figure 32. This would be very obvious if one were to plot the absolute values of the stresq the actual values are of course negative because the stress is compressive. As noted above, this causes the rate constant for the silylation reaction to decrease with time, and this (together with the increased diffu&ity of P) leads to a readion front that is not as sharp as it would otherwisebe. Finally, Figures 35 and 36 show (as solid curves) the overall layer thickness 1 as a function of time. The time span in Figure 35 is the same as in Figures 29-34. Since the value of 1 at T = 10 is still less than 1.4, the process is obviously very far from completion (1 = 2.0), in accordance with Figure 32. Figure 35 also gives results for a difIerent value of Poisson’s ratio, namely v = 0.3, which is appropriate to a wide variety of elastic solids. The volume increasein this case is considerably smalleq this is expected, because the compressivestresseswill have a larger effect on the volume the more v deviates horn the value of 0.5 for an incompressiblesolid. The behavior of l(~) over the entire course of the process is shown in Figure 36. There are four identifiable although indistinct regimes: a short induction period (more obvious in Figure 35) due to the preliminary formation of U; a period of rapid, unconstrainedexpansio~ a period of slowergrowth after the leading edge of the expansion front reachesthe impermeableboundary; and a iinal stage of no growth after the reactions are complete. The uMrnate vahe of 1 at T = cm can be obtained from Equation (185) by reakzing that the layer thicknessis inverselyproportional to the total concentration of Q, U, and E at any time. This gives 1= =I+(l+v)(a-l) 3(1 – v) (195) For v = ~ the result is simply lW = a, as expected. For v = 0.3 and a = 2 one finds lm = 1.619, which agrees well with the asymptoticvalue in Figure 36. Many more piwuneter studies could be undertaken,but it should be rememberedthat the resultsme largely fictitious due to the assumption of purely (and linearly) elastic behavior. Therefore, as noted earlier, this model should be used primarily to check the functioning of the multidimensionalcodes. 63164 (a) Cross-section prior to TSI (c) Post exposure bake ,. (e) Plasma descum (f) Plasma etch Figure 1: TSI processing steps. 65 unsilylated resist \ \ silicon wafer substrate ~ Figure 2: Cross-sectionalplasma stain of silylationprior to etch . . /w va Figure 3: Arbitrary moving control volume. 67 (t) arb ●❐ Vp(t) Figure 4 A two-dimensionalquadrilateralcontrol volume centered about point P. 68 I m r --.-.90 ------ I k-------------- A9 gas, g ------- -- ~Ay ~ -...---k--J surface, s resist, r I ■ mA r Figure 5: Control volume and fluxes at resist/gas interface. 69 .— B“= .-. -. W ,- .— , / / / ,— 4 /zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA /’%. / ,.% ...’ ,-x.~ ‘ ,..-.,..’ ,-.. ”” ,-, / A... .0 .<; .,’’-’” -’ ‘--.,-’. ..-, ,“<.- / ~------ .--=- ~/ —--- ~ - -—’-’- / Figure 6: The Eyring dashpot element (upper figure) and the Langevin spring element (lower figure). .. . .... . 12 1 I .. I . zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDC I I : . 10 . . m * ■ . 8 . . m m m 6 . . 4 . . 2 . , , m 0 0 5 10 15 20 25 Stress ( x 107dynes/cm2) Figure 7 Dashpot response. 71 30 35 Stress ( x 107dynes/cm2) o N o 8 ~ i5 ~ o * ● O . M * o b . * f!” @l .0 al . * * o @ . , A , Y 3 4 1 x Figure 9: Four node qu~aterd elemmt in giobd coortite 73 systa, -1 1 F@re 10: Four node quadrilateralelementin local isoparametnc coordinate system. 74 3 4 fxmp4’ymp4) -m mm- ----- =---- 1 (%pl’ympl 2 , Figure 11: A four node quadrilateral element showing sub control volumes and integration points for the CVFEM. 75 5X108 4xlo8 3X108 I 2)tlo8 lxlo8 o 0 0.1 0.05 ‘TrueStrain Figure 12: Response in uniaxkd extension. 76 0.15 0.2 5xlo8 ~lrwo. -“ 111,,,,, ‘1111, , / :* j 4xio8 -- “’’1 s11,,,,,,,, zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQ 1111111 * 1 1 ,1 1 ,1 1 ,,8 1 , ‘* ,14,181,118111 N E ‘~-~m- t) 2 g 3xlo8 lDl=O.1~ %’ Note: ‘D’=~ IX108 ID Ik a o 0 proportional to strain a B 0.05 0.1 True Strain F@e 13: Rate dependence of flow strength in uniaxkd extension. 77 rate 1 0.15 5X108 I “~ “% h = 0.9 E9 dynes/cm2 * 9 4X108 .. N ,. I I \ - +* “,” . . . . . . . . . . , ,, .,,,,,,,,,, h = 9.0 E9 dynes/cm2 =,% r ---*..-,,, ,,, ,,, ,,, , #,, .,,.,,.,,. -,,, ,, .,,.,,,, ,,,.-,,, ,,, ,m ,,7 / h = 90. E9 dynes/cm2 - 1X108 1 o 0 0.05 I 0.1 I 0.15 True Strain Figure 14: Softening response as a function of modulus in uniaxial extension. 0.2 . 6xI08 . I 1 I . 5X108 4X108 ,11111,,,,,,8,,+, ,, , , 3X108 r S~~=62 E7 dynes/cm 2 % * 2X108 m S,.= 50 E7 dynes/cm 2 * s m ■ . 1X108 . I o 0 0.05 I I 0.1 ‘0.15” 0.2 True Strain Figure 15: Steady state flow strength response in uniaxial extension. , 6X108 F 1 I I 5X108 m 1 n ~ 1 1 m 1 I 1 I Material Point Simulation m 1 1 1 1 2D Model ~ ■ a 1 1X108 m # 1 1 o o 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 True Strain Figure 16: Large strain uniaxial extension simulation: Comparison of material point simulation and 2D model. 80 0I I i 1 8 J 1 1 I 1 I , # 1 I \ -1)(108 0 I I 1 I 1 I I 1 I 1 1 , I 8 I I { YY -2)(108 -3.108 -i -4.108 -5X108 True Strain Figure 17: Response in free surface swelling. 8.8x108 8.6X108 . 8.4x1 08 . 8.2x108 . \ \ 8X108 , 7.8x108 t o . . .. u. . . . . . . . . . . . . . . . . . . . 0.05 0.1 0.15 True Strain Figure 18: Strain softening response in free surface swelliig. 82 0.2 0.25 0.35 0.3 0.25 u) 0.2 0.15 0.15 0.2 0.25 0.3 0.35 0.4 Poisson’s Ratio Figure 19: Free surface expansion for ideahzed elastic response. 83 0.45 0.5 0.3 1 s 1 b 1 m n 8 1 0.25 m Material Point Simulation 1 1 t 1 m 0.2 1 ■ m m \ 1 1 c “q 1 0.15 # m m 1 G 1 n # n 0.1 1 D m I I n n n 0.05 1 m m E 1 n I 1 1 n 1 0 0 2X1O‘6 4X1o ‘6 6X10-6 8x1O‘6 IXIO-5 Time -s Figure 20: Large strain swellingsimulation: Comparison of material point simulation and 2D model. Free end boundary conditions specified species concentrations (J=O / 100 elements Symmetry boundary conditions - mass - fluxes = O - Yk y% \ end boundary conditions mass fluxes = O Fixed Cl=v=o F@re 21: Geometry for 1-D solution comparisons. 85 2-D Model (~=.5), \ 2.10-5 i # 1.8x10-5 / b4 4’7 2-D Model 4 Elastic/plastic 4 “1 2-D ~u 1.6x10-5 (V=.3) Model(V=.3) - ❑ -mo u 9-m ■ -m Model (V=.3) 1 .4X1 0-5 ■ 1.2 X10-5 r IXIO-5 o 20 40 60 80 100 Time -s Figure 22: 1-D and 2-D model predictionsfor layer thicknessas a function of time. 120 ●ooool-D 1.2 — J Model 2-D Model t=2.5 s ~t=60 S 0.6 0.8 =2.5 S 0 o 0.2 0.4 Normalized 0.2 o 1 0.4 Normalized length (z) 0.6 length(z) 0.8 1 (b) Distribution of unsilylated resist (a) Distribution of silylating agent o * I I.... — Model 2-DModel 1-D ● -2X109 -4X109 “E < ~ -6x109 D ; ~ -8x109 3 . . . . . 1-D Model - t=2.5 s — 2-D Model - -Ixlo’” 0 0.2 0.4 0.6 0.8 -1.2X10’C o 1 Normalizedlength(z) I 0.2 I 0.4 I 0.6 t 0.8 Normalizedlength(z) (c) Distribution of silylated resist (d) Lateral stress distribution F@re 23: Distributions of species concentrations and stress at three diiTerenttimes for v = 0.5. 87 1 NtMe: dimensions in cm t--–+, 2“4’’0-5 8)(10-6 ~– ~1 1 Free sutiace -f- CrOsslin~w! u ncrosslin~ed x Fixed surface F@re 24: Computational mesh for the nominal 2D silylation simulation. 88 HI seconds E!El t=l second k3 seconds Mass Concentration Species A in @cm3 1! — — t=5seconds 2. lzE+OO 2.13 EIE+EVZi 1. 88E+Ela 1. 75E+EWI 1. 63E+EEi 1. 5BE+BD 1. 37E+EIE! 1. 25E+12EI 1. 13E+~B 1. WZ!E+EIEI 8.75E-01 7.5DE-~1 6.25E-E)I 5.WIE-EI1 3.75E-EIl 2.512E-BI 1.25E-EII 0.00E+OO 1 t=10seconds ii mkmdiia t=30seconds t== seconds F@ure25: 2Dmatetid sweMngtith supertiposed concentration pro~esof thestiylation agent (species A) as a function of time. 89 ~,.—_= . Species B concentration glcms ti.;% -—-. .—.. b Species C concentration glcms .2.mz-01 1. 28E+EE i.e8E-el 1.2EE+E0 1.13E+E0 1.85E+EB 1.7sE-al 1.63Z-81 1.50E-01 1,., 1.3sE-al 1.Z5E-01 1.12E-’31 1.eEE-al e.75E-’JZ 7.sOE-’a2 6.25E-E2 5.eeE-E2 3.75E-EZ 2.5EE-132 1..?5C-E2 E.EEE+OO -1.25E-02 Species D concentration glcms I 9.75E-E1 9.eBE-el S.25E-L41 7.5aE-El 6.75E-01 6.8EE-01 1 5.25E-Eil 4.SEE-U 3.75E-01 3.8EE-01 2.25E-01 1.50E-91 7.5EE-E2 ~.EEE+W Species E concentration glcms 3.2EE-E2 3.EOE-82 2.eEE-E2 2.6’3E-02 2.4BE-W 2.2EE-E2 2.EaE-02 1.eEE-a2 1.6EE-~2 1.4EE-i32 1.2EZ-B2 1.‘3aFz-e2 a.EEE-e3 6.8~E-03 4.EEE-E3 2.EEE-E3 B.aaz+m -.2. EaE-03 Figure 26: 2D concentrationprofilesfor species B, C, D, and Eat 60 seconds. ..----------- Sn .-___ -- .... ..m.-- . —— . . . :.. . dynes/cm2 ‘: *._... — -- .--— a “ ‘“ J ,. ,. . . . —-. .. .. . - Syy dynes/cm 2 4 —.. ,/ -,—. .1 . ., @ >. .. Sw dyneslcm2 / -... \ ‘@’ . .. --—-. Figure 27: . . ‘t--- -- .: ? ,, ,, 7’ 2D stress dktributions at 60 seconds. 91 CWx3slinkeciiekments UuwwWnkeci elementsf Figure 28: 2D deformed mesh at 60 seconds. 92 7 1. c .-o % 0. L Ea) 0 c 0 0. 0 (n m a) z 0 .- 0. tn c .E— n 0. 0. 0.0 0.2 0.4 0.6 Dimensionless position 0.8 1.0 z Figure 29: Concentration profilesfor silylatingagent P in one-dimensionalmodel for O= 0.01, P = 1, T. = 1, a=2,7=l, v = 0.5, and A = 1. 93 1. 0 0. 8 0. 6 0.zyxwvutsrqponmlkjihgfedcbaZYXWVUTSRQPONMLKJIHGFEDCBA 4 time o.2- 0.o0.0 1 0.2 I 0.4 1 0;6 Dimensionless position i 0:8 f 1:0 z Figure 30: Concentration profilesfor unreactedsolid Q in one-dimensionalmodel for 0 = 0.01, ~ = 1, T, = 1, Cr=z,’y= l,v=0.5, and A=l. 94 0.3 k= 0.5 1.0 0.2 0.1 0.0 0.0 0.2 0.4 0.6 Dimensionless position 0.8 1.0 z Figure 31: Concentration proiilesfor reacted but unexpanded solid U in one-dimensionalmodel for O= 0.01, ~=1, ~,=1, cx=2, -y=l, v= 0.5, and ~=1. 95 : a) z 0 .— CA c a) E .— n 0.0 0.2 0.4 Dimensionless 0.6 position 0.8 1.0 z Figure 32: Concentrationprofilesfor reacted and expandedsolid E in one-dimensionalmodelfor O= 0.01, P=l, ~r=l, a=2, y=l, v=0.5, and A=l. 96 0.05 -e- 2.0 3.0 , 1.0 0.04 % .- 0 0 : 0.03 % a) F .-0 0.02 (I3 c .-? n 0.01 0.00 0.0 0.2 0.4 Dimensionless 0.6 position 0.8 1.0 z Figure 33: Velocity profiles in one-d~ensions,l model for O= 0.01, ~ = 1, r, =1, a =2, ~ = 1, v = 0.5, and A=L 97 o.,0 -o. 1 -o. 2 -o. 3 time -o. 4 0.0 0.2 0.4 Dimensionless 0.6 position 0.8 1.0 z Figure 34: Stress profiles in one-dimensionalmodel for O= 0.01, ~ = 1, I-, = 1, a = 2, ~ = 1, v = 0.5, and A=l. 98 . ... 1.4 1 1.3V=O.5 1.2u) u) a) T 0 .— rn c .-; n 1.1 - 1,0 ‘~ 0 I 2 6 4 Dimensionless time Figure 35: Progress of swelling in one-dimensional model for O = 0.01, @=l, A=l. 99 10 8 ~r=l, a=2,’y=l, and 2.0 1.8 V=O.5 .-0---===--= =.==. =-=..... ...................... 1.6 1.4 1.2 1.0 o 20 40 Dimensionless 60 80 100 time Figure 36: Overall progress of swellingin one-dimensionalmodel for 6 = 0.01, /3= 1, 7. = 1, a =2, ~ = 1, and~=l. 100 References [1] C. Henderson, D. Wheeler, T. Pollagi, G. Cardinale, D. O’Connell, A. F~her, and J. Goldsmith. “Top Surface Imaging for Extreme Ultraviolet Lithography”,. Y. Vladimirski, editor, Emerging Lithographic TechnologiesI.. Proc. SPIE 3331, 1998. [2] W.S. Han, J.H. Lee, H.Y. Kang, C.G. Park, Y.B. Koh, and M.Y. Lee. “Optimization of Deep UV Positive Tone Top SurfaceImaging Process”,. Microelectronic Engineering, 23:255-258,1994. [3] D.C. La Tulipe, A.T.S. Pomerene, J.P. Siions, and D.E. Seeger. “Positive Mode Silylation Process Characterization”,. Microeimtronic Engineering, 17265-268, 1991. [4] M.W. Horn, B.E. Maxwell,R.B. Goodman, R.R. Kunz, and L.M. Eriksen. “Plasma-depositedSilylation Resii for 193 nm Lithography”,. J. Vat. Sci. Technol. B, 14(6):4207-4211, 1996. [5] T. Itani, H. Yohino, S. Hashimoto, M. Yamana, N. Samoto, and K. Kasama. “A Study of Acid Diilusion in Chemically AmplMed Deep UltravioletResist”,. J. Vat. $ci: l%chnol.B, 14(6):4226-4228,1996. [6] N. Glezos, G.P. Patsis,I. Raptis, and P. Argitis. ‘Application of a Reaction-Dfilon Model for Negative Chemically knplified Resiits to DetermineElectron-Beam Proximity Correction Parameters”,. J. Vat. Sci. Technol. B, 14(6):4252-4256, 1996. [fl C.S. Whelan, D.M. Tanenbaum, D.C. La Thlipe, M. Isaacson, and H.G. Craighea.d. “Low Energy Electron Beam Top Surface Image Processing Using Chemically Amplified AXT Wit”,. J. Vat. Sci. Technol. B, 15(6):2555-2560, 1997. [8] M.A. Zuniga and A. R. Neureuther. ‘Diagnostics of Patterning Mechanismsin Chemically Ampliiied Resiits from Bake Dependencies of Images”,. J. Vat. Sci. Technol. B, 14(6):4221-4225, 1996. [9] L.E. Ocola and F. Cerrina. “Latent Image Characterizationof Postexposure Bake Process in Chemically Ampliied Resists”,. J. Vat. Sci. Technol. B, 15(6):2545-2549, 1997. [10] M. Weiss and M. Goethals. “A Novel Approach to Simulation of the Silylation Bake in the DESIRE Process”,. Microelectronic Engineering, 30:313-316,1996. [11] B.E. Deal and A.S. Grove. “General Relationship for the Thermal Oxidation of Siicon”,. ~hgS., J. Appl. 36(12):3770, 1965. [12] L. Bauch, U. Jagdhold, J. Bauer, B. Dietrich, and W. Hoppner. ‘Positive Mode Siiylation Process Characterization”,. MicroelectronicEngineering, 17265-268,1991. [13] M.A. Hartney, M. Rothschild, R-R. Kunz, D.J. Ehrlich, and D.C- Shaver. “PositiveTone Siiylation Processes at 193 nm”,. Microelectronic Enginee%ng, 1351–56, 1991. [14] P.J. Paniez, O.P. Joubert, J.C. Oberlin M.J. Pens, T.G. Vachette, A.P. Weill, J.H. Pelletier, and C. FIori. “The Role of Self Dfion in the Dry Development and Plasma Durability of Polymers”. Microelectronic Engineering, 13:57-60,1991. [15] C. Pierrat. “New Model of Polymer Silylatioxx Application to Lithography”,. J. Vat. Sti. Technd B, 10(6):2581-2588, 1992. [16] W.S. Winters and W.E. Mason. “The Computer Code PASTA and Its Application to Coupled Thermal- 101 Mechanical Problems”,. S.N. Atluri and G. Yagawa, editors, The Proceedings of the International Conference on Computational En~”neering Science, pp. 59.vii.2-59.viL4, Springer-Verlag, New York, NY. 1988. [17] W.S. Winters and W.E. Mason. ‘PASTA2D-Applications”,. Technical Report SAND89-8217, San&a National Laboratories, Livermore, CA, 1989. [18] Y.L Dimitrienko. “Effect of Ftite Deformationson Internal Heat-Mass Transferin Elastomer Ablating Materials”,. Int. J. Heat Mass fiansfer, 40(3):699-709, 1997. [19] M.A.Zuniga and A. R,.Neureuther. ‘stress I)ependent Sfiylation Model and Tw~Dimensional Profile Simulation”,. J. Vat. Sci. T&not. B, 15(6):2565-2569, 1997. [20] D. Wheeler, E. Sciwwrer,G. Kubiak, R. Hutton, S. Stein, R. Cirelli, F. Baiocchi, M. Cheng, C. Boyce, ad G. Taylor. “SilylatingReagents with High Siicon Contentsfor Dry-Developed Positiv&Tone Resists for Extreme-W (13.5 nm) and Deep-UV (248 nrn) Microlithography”,. R. L. Clough and W. Shalaby, editors, Irradiation of Polymers: Fundamentalsand Technologi~l Applications; A G’SSymposiumSeries No. 620, pp. 399-415. The American Chemical Society, 1996. [21] A. Wong. PhD. Thesis. The Universityof California, Berkeley, CA, 1994. [22] S. Whitaker. FundamentalPrinciples of Heat T#-ansfer.R. E. Krieger PublishingCo., MaLabor, Florida, 1983. [23] M. A. Hartney. ‘Modeling of Positiv~tone Silylation Processes for 193-nm Lithography”,. J. Vat. Sci. Technol. B, 11(3):681-687, 1993. [24] N. L. Thomas and A. H. Wmdle. “A Theory of Case II Dif%@on”,. Polymer, 23:529-542, April 1982. [25] J. Crank and G. S. Park. D@sion in Poiymers. Academic Press, New York, 1968. [26] C. Gostoli and G. C. Sarti. “DifhAon and Localized SwellingResistancesin GlassyPolymers”,, Polymer Engineering and Science, 22(16):52%542, 1982. [27] J. Crank. The Mathematics of Di..ion, 2nd Edition. Clasendon Press, Oxford, 1975. [28] S. Middleman. Introduction to Mass and Heat T?ansfer. Wiley, New York, 1998. [29] R.N. Haward and G. Thackray. “The Use of a MathematicalModel to Describe IsothermalStress-Strain Curves in Glassy Polymers”,. Proc. Roy. Sot. A., 302:453-472,1968. [30] M.C. Boyce, D.M. Parks, and A.S. Argon. ‘Large Inelastic Deformation of Glassy polymers. Part I Rate Dependent Constitutive Model”,. Mechanics of Materials, 7:15-33, 1988. [31] A. Hoger. “The Material Time Derivativeof Logarithmic Strain”,. Int. J. solids Structures, 22(9):1019- 1032,1986. [32] D.J. Bammann and G.C. Johnson. “On the Kinematics of Finite Deformation Plasticity”,. Acts Mwhanicu, 70:1-13,1987. [33] D.J. Bammann and E.C. Aifantis. ‘A Model for Finite Deformation Plasticity”,. Acts Mechanicu, 69:91-11711987. [34] K.H. IIuebner and B.A. Thornton. The Finite Element Method for Engineers. John Wiley, New York, 102 1982. [35] P. N. Brown, G. D. B~e, and A. C. Hindmarsh. “VODE: A Variable CoefficientODE Solver”,. SIAM J. Sci. Stat. Comput., 10:1038-1051, 1989;. also Lawrence Llvermore National Laboratory (LLNL) Report UCRL-98412, June 1988. [36] G. E. Schneider and M. J. Raw. “Control Volume Ftit&Element Method for Heat Transferand Fluid Flow using Colocated Variables-1. Computational Procedure”,. NumericalHeat tinsfer, 11:363-390, 1987. [37] B. F. Blackwell and R. E. Hogan. “Numerical Solution of Axiaymmetric Heat Conduction Problems using Ftite Control Volume Technique”,. Journal of Thermophysics and Heut Transfer,7(3):462471, 1993. [38] R. E. Hogan, B. F. Blackwell,and R. J. Cochran. “Application of Moving Grid Control Volume Finite Element Method to Ablation Problems”. Journal of Themtophysics and Heat lhmsfer, 10(2):312-319, 1996. [39] J.O. Hallquist. “User’s Manual for DYNA2D - An Explicit Two-Dimensional Hydrodynamic Ftite Element Code with Interactive Rezoning”,. Technical Report UCID-18756, Lawrence Llvermore National Laboratory, Livermore, CA, 1980. [40] L.M. Taylor and D. Flanagan. “PRONT02D, A Two-DimensionalTransientSolid Dyntics Program”,. Technical Report SAND860594, Sandia National Laboratories, Albuquerque, NM, 1987. [41] J.O. Hallquist. “NIKE2D: An Implicit, Ftite-Element Code for Analyzing the Static and Dynamic Response of Two-Dimensional Solids”,. Technical Report UCRL52678, Lawrence Llvermore National Laboratory, Llvermore, CA, 1979. [42] G. Maenchen and S. Sack. ‘“j?heTENSOR Code”,. Methods in Comp. Physics, 3:181-210,1964. [43] D.P. Flanagan and T. Belytschko. ‘A Uniform Strain Hexahedron and Quadrilateralwith Orthogonal Hourclass Control”,. Int. Journalfor Numerical Methods in Engineering, 17679-706,1981. [44] T. Belytschkoand J.S. Ong. ‘Hourglass Control in Linearand NonlinearProblems”,. Computer Methods in Applied Mechanics and Engineering,43:251–276, 1984. [45] K.J. Perano and V.N. Kaliakin. “INTERP- A Fortran Callable Free Format Data Interpretation Subroutine System”,. Technical Report SAND87-8244, Sandia National Laboratories, Lhwrmore, CA, 1989. [46] K.J. Perano. “COMPRO - A SubroutineSystem for Syntactical Analysis”,. TechnicalReport SAND898441, Sandia National Laboratories, Livermore, CA, 1989. [47] P.E. Nielan, K.J. Perano, and W.E. Mason. “ANTIPASTO: An Interactive Mesh Generator and Preprocesmr for Two-Dimensional Analysis Programs”,. Technical RepOrt SAND90-8203, Sandia National Laboratories, Livermore, CA, 1990. [48] R.El. Bird, W.E. Stewart, and E.N. LightfixX. TkansportPhenomena. Wley, New York, 1960. [49]L.D. Landau and E.M. Ltitz. 17WOW oj Ekwticity.PergamonPress, Oxford,1970. 103 [50] .J.R.Welty, C.E. Wkks, and R.E. Wilson. Fundamentals of Momentum, Heat and Mass Thansfer. Wiley, New York, 1969. [51] L.R. Petzold. “A Description of DASSL: A Differential-AlgebraicSystem Solver”,. Technical Report SAND82-8637, Sandia National Laboratories, Livermore, CA, 1982. 104 .. ,.- .— —-.--- .- . . . . . Distribution EXTERNAL DISTRIBUTION: Prof. D. S. Dandy Colorado State University Dept. Agriculture and Chem Eng. Fort Collins, CO 80523 Prof. M. C. Boyce Massachusetts Institute of Technology Dept. of Mechanical Engineering Cambridge, MA 02139 L. Chuzhoy Caterpillar Inc. Technical Center, Bldg. E Division 854 P.O. Box 1875 Peori~ IL 61656-1875 J. Cobb EWLLCiMotorola Assignee Sandia National Laboratories P.O. Box 969 MS 9911 Livermore, CA 94551-0969 E. Croffle 211-19 Cory Hall #1772 EECS Dept. University of California Berkeley, CA 94720 Prof. H. A. Dwyer Dept of Aero. and Mech. Eng. University of California Davis, CA 95616 105 .- Prof. R.W. Dutton CIS Extension, Room 333 Integrated Circuits Laboratory Stadord University Stanford, CA 94305-4075 G. Feit SEMATECH 2706 Montopolis Drive Austin, TX 78741-6299 Prof. D. L. Flamrn Dept. EECS 187M Cory Hall lhiversi~ of California Berkeley,CA 94596. Prof. R. Greif Dept. of Mech. Eng. University of California Berkeley, CA 94720 Prof. R. J. Kee Engineering Division Colorado School of Mines Golden, CO 80401-1887 Ellen Meeks (RD) Reaction Design 6500 Dublin Blvd., Suite 214 Dublin, CA 94568 Prof. A. R. Neureuther Dept. EECS University of California Berkeley, CA 94720 Uzo Okoroanyanwu One AMD Place P.O. Box 3453 MS 78 Sunnyvale, CA 94088-3453 106 Prof.W. G. Oldham Dept.EECS University of California Berkeley, CA 94720 Dr. S. C. Palrnateer Room MEL 204B 244 Wood St. Lexington, MA 02173 Veena Rao Intel SCI-03 3065 Bowers Ave. Santa Clara, CA 95052 G. Taylor Shipley 455 Forest St. Marlborough, MA 01752 C. Grant Willson Chemical Engineering Campus Mail Code: C0400 Universi@ of Texas Austin, TX 78712 Marco Zuniga Clo Prof. A. R. Neureuther Dept. EECS University of California Berkeley, CA 94720 107 INTERNAL DISTRIBUTION: MS 1405 MS 9001 1126 1126 1711 1812 8000 MS 9214 MS 9420 8980 8200 MS 9405 MS 9409 MS 9409 MS 9409 8230 8230 8250 8250 8270 8270 MS 9054 8300 MS 0601 MS 0601 MS 0603 MS9405 MS 9409 MS 9041 MS 9042 MS 9051 MS 9042 MS 9042 MS 9051 MS 9051 MS 9042 MS 9042 MS 9402 MS 9051 MS 9042 8345 8345 8345 8345 8345 8345 8345 8345 8345 8345 8345 8345 M. E. Coltrin J. Y. Tsao P. Esherick D. R. Wheeler T. O. Hunter Attn: J. B. Wright, 2200 M. E. John, 8100 R. C. Wayne, 8400 L. A. Hiles, 8800 D. R. Henson, 8900 W. E. Mason A. L. West Attn: L. N. Tallerico, 8204 A. J. West, 8240 C. T. Oien, 8260 J. M. Hruby C. C. Henderson G. D. Kubiak G. F. Cardinale R. H. Stulen A. Ray-Chaudhuri W. J. Mclean Attn: L. A. Rahn,8351 F. P. Tully,8353 W. Bauer, 8358 D. R. Hardesty, 8361 R. W. Carling, 8362 R. J. Gallagher, 8366 J. S. Binkley G. H. Evans (10) J. F. @Car s. Grifflths W. G. Houf R. S. Larson A. E. Lutz C. D. Moen R. H. Nilson D. Rader F. M. Rupley P. A. Spence 108 — MS 9051 MS 9042 MS 9052 MS 9405 8345 8345 8361 8700 A. Ting W. S. Winters (10) M. Allendorf M. T. Dyer Attn: M. I. Baskes, 8712 J. C. F. Wang, 8713 G. J. Thomas, 8715 K. L. Wilson,8716 S. M. Foiles,8717 W. A. Kawahara,8746 MS 9402 MS 9042 MS 9405 8701 8742 8743 C. M. Hartwig E. Chen P. E. Nielan Attn: D. J. Bammann,8743 M. Horstemeyer, 8743 M. L. Chiesa, 8743 J. J. Dike, 8743 A. R. Ortega, 8743 MS 9405 MS 9405 MS 9405 MS9011 MS 0841 8743 8743 8743 8980 9100 L. A. Bertram J. F. Lathrop V. C. Prantil (10) K. Perano P. L. Hommert Attn: A.Ratzel,9112 R. Griffith, 9114 W. Rutledge,9115 C. Peterson,9116 T.Baca,9119 I MS 0826 MS 0826 MS 0826 MS 0826 MS 0443 9111 9111 9111 9113 9117 MS 0443 MS 0443 MS 0443 MS 0443 MS 0443 9117 9117 9117 9117 9117 W. L. Hermina D. K. Gartling R. Shunk S. Kempka H. S. Morgan Attn: J. B. Aidun, 9117 S.W.Attaway,9117 S.N.Burchett,9117 A.F.Fossum,9117 M. L. Blanford S. W. Key M. K. Neilsen V. L. Porter G. W. Wellman 109 3 1 1 MS 9018 MS 0899 MS 9021 1 1 MS 9021 MS 1380 Central Technical Files, 8940-2 Technical Library, 4916 Technical Communications Department, 8815/ Technical Library, MS 0899,4916 Technical Communications Department, 8815 For DOE/OSTI Technology Transfer, 4212 110