Abstract:
The invention relates to the simulation method and apparatus used in plasma modeling. It includes a method to transform transient formulations of the phenomenological plasma model into a quasi-stochastic spatial formulation. Specifically, the invention aids in decreasing computational time for the modeling of plasma in a plasma processing system, particularly those involving two different time-based parameters. The invention is particularly described in connection with plasma simulations used for the optimization dual-frequency capacitively-coupled plasma etching systems.

Description:
CROSS-REFERENCE TO RELATED APPLICATIONS 
       [0001]    This application claims the benefit of U.S. Provisional Application Serial. No. 61/094,719, filed on Sep. 5, 2008, the disclosure of which is hereby incorporated by reference herein in its entirety. 
     
    
     FIELD OF THE INVENTION 
       [0002]    The invention relates to the modeling of plasma distribution properties in the design and optimization of plasma processing systems. More particularly, the invention relates to a system and a method for modeling the distribution within the chamber of a plasma processing system, of plasma density, plasma species, or other plasma properties. The invention is particularly useful in the design and optimization of dual-frequency CCP etching systems, IPVD deposition systems, and other plasma processing or plasma enhanced processing systems. 
       BACKGROUND 
       [0003]    Advanced plasma tools have been using very high frequency power sources or combination radiofrequency (RF) and direct current (DC) power supplies to advance greater control over the plasma formation process and performance. Recently, for example, capacitively coupled plasma (CCP) reactors have employed dual frequency (DF) power delivery to control the energy and density of the ions and radicals at the surface of a wafer being processed. In such reactors, the effectiveness and the quality of the process being performed depends on the distributions of this energy and of the various particle species in the processing chamber, particularly near the wafer surface. Inductively coupled plasma (ICP) reactors, such as ionized physical vapor deposition (iPVD) reactors, employ combinations of gas and coating material species plus energy introduced from a number of different sources within the geometry of a processing chamber. In these systems, the material and energy distributions that affect the process being performed are reached after initial parameters are established, transient behavior has subsided, and the system has stabilized. However, the physics of this process has not been fully understood and has been the subject of continued theoretical and experimental investigations. 
         [0004]    Plasma modeling and numerical simulations play in an important role in understanding plasma behavior in these systems. Modeling and simulation methods can play important roles in the development and optimization of plasma equipment and related plasma processes. 
         [0005]    While, in the last decade, computing systems have greatly increased capabilities for numerically modeling complex and robust plasma-chemical systems like etching systems (for example central processing unit (CPU) speed, multi-core architecture, memory resources, and large size data storage capacities), the computational time required for a full numerical simulation is often on the scale of many hours, and in some specific complex cases, can be several days. Most computing systems and software are not yet ready and/or capable for parallel computing for a variety and specific tasks. 
         [0006]    For example, even with the simplest plasma model, e.g. inert argon gas with a basic geometry etching system, a computer simulation of the plasma density distribution can take several hours (2-6 hours) to converge if the model includes mutually coupled mass transport equations, heat and energy balance equations, and electromagnetic (EM) wave propagation equations of 100 MHz (i.e. RF cycle on a 10 ns time scale). To generate a full transient solution on a system excited by a lower frequency EM wave would require yet longer CPU times. 
         [0007]    Computational time becomes even longer when more complex chemistry (number of reactant species and/or number of participating reactions) and higher resolution reactor geometry are modeled. Computing time increases yet again if a transient solution is desired, when a pulsed process is investigated by simulation, or when a process is operating on multiple frequency scales, such as in the DF-CCP etching system. Multiple frequency scales or other multiple time based factors increase computational time because electron and ion transport properties within the plasma must be solved at both time scales. For example, the first frequency of the DF-CCP etching system can be 100 MHz (corresponding to a 10 ns time scale); and the second frequency can be 1-2 MHz (corresponding to a 1 μs time scale). Thus, the DF-CCP etching system must be solved for both the 10 ns and the 500 ns (or 1 μs) time scales. The computational time is related to the product of the number of iterations that each solution requires. 
         [0008]    If the interactions between the ion flux and the processed wafer are considered, which can be over the time scale of several thousand RF cycles, the number of iterations and computational times can increase drastically. 
         [0009]    Thus, there is a need to develop fast methods that would reduce computational real time in complex chemistry and geometry plasma systems. The fast algorithm and/or method could apply to solving numerically by finite elements modeling (FEM) methods, multiple frequency applications, pulsed plasma processes, sequential processing, contamination studies, feature profile evolutions, or developing industrial applications for etching or deposition technologies in plasmas. 
         [0010]    Further, there is a need to develop and invent new, fast algorithms for technological plasma simulations and modeling. The simulation procedure should be practical and efficient for aiding in process and/or hardware development, for serving as an operational and predictive flexible tool, and for generating outputs as close as possible to the development needs and workflow in real time. Thousands of wafers are processed by the technological process on a daily basis, thus a large amount of experimental data related to hardware and software performance is readily available. Simulation procedures should be competitive in the rate of generating results to provide the ability to correlate, correct, or change the technological process in real time. This can represent large material and technical resource savings. 
       SUMMARY OF THE INVENTION 
       [0011]    According to principles of the present invention, a computer method is provided for simulating plasma density distribution within a plasma reactor that is responsive to at least two variables of different time bases. The computer method includes executing computer program codes according to a process of approximating transient sheath behavior of the plasma occurring in the plasma reactor in response to a variable having a first time basis and performing, on a spatial character that is based upon the probability function, a time-domain simulation of the transient plasma sheath behavior occurring in the plasma reactor in response to a variable having a second time basis. The behavior is approximated by defining a model domain that is characteristic of the plasma reactor, constructing a probability function for a time dependent variable within the model domain, and transferring the time dependent variable into the spatial character. In certain embodiments, the performing is carried out to arrive at a steady state plasma density distribution. 
         [0012]    According to an exemplary embodiment of the invention, a computer method is provided for simulating plasma density distribution within a capacitively coupled plasma reactor having a dual frequency power source coupled to the plasma, where transient sheath behavior of the plasma occurring in the plasma reactor at a first frequency of the dual frequency power source is approximated and a time-domain simulation of the transient plasma sheath behavior occurring in the plasma reactor is performed on the spatial character at a second frequency of the dual frequency power source. 
         [0013]    In certain embodiments of the invention, the constructing includes random time sampling of the time dependent variable and constructing a probability function of a predetermined plasma parameter response to an incremental unit of a time dependent variable within the model domain. In some embodiments of the invention, the defining includes establishing at least one initial value and at least one boundary condition for the model domain. 
         [0014]    A plasma reactor optimized in accordance with the simulation performed by the computer method of the invention is also provided. For example, a dual-frequency plasma etching reactor optimized in accordance with the simulation performed by the computer method of the invention is provided. 
         [0015]    A computer program product is also provided having computer readable program codes embodied in the medium that when executed causes a computer to execute the method of the invention. A simulation apparatus is also provided having a calculating unit programmed with the computer program product. 
     
    
     
       BRIEF DESCRIPTION OF THE FIGURES 
         [0016]      FIG. 1  is a diagram of a computer system environment for the present invention. 
           [0017]      FIG. 2  is a cross-sectional diagram of a DF-CPP etching system to be modeled in accordance with the description of an exemplary embodiment of the invention. 
           [0018]      FIG. 3  is a flow chart and information flow diagram of a simulation or modeling process for solving the transient plasma density distribution in a plasma processing system such as, for example, the DF-CCP etching system of  FIG. 2 . 
           [0019]      FIGS. 4A and 4B  are flow charts of a simulation and modeling process using a stochastic transformation according to one embodiment of the present invention. 
           [0020]      FIG. 5  is a schematic illustration of the transient character of the points within the simulation domain. 
           [0021]      FIG. 6  is a schematic illustration of the characteristic values and the arrangement of the characteristic values into arrays, B n . 
           [0022]      FIG. 7  is a schematic illustration of the random sampling of Ω for each sub-domain within the stochastic simulation domain. 
       
    
    
     DETAILED DESCRIPTION 
       [0023]    Turning to the drawings, wherein like numbers denote like parts throughout the several views,  FIG. 1  illustrates a hardware and software environment for a computing system  10  that may include an integrated circuit device (hereinafter “ICD”) consistent with embodiments of the invention. The computing system  10 , for purposes of this invention, may represent any type of computer, computer system, computing system, server, disk array, or programmable device such as multi-user computers, single-user computers, handheld devices, networked devices, etc. The computing system  10  may be implemented using one or more networked computers, e.g., in a cluster or other distributed computing system. The computing system  10  will be referred to as “computer” for brevity sake, although it should be appreciated that the term “computing system” may also include other suitable programmable electronic devices consistent with embodiments of the invention. 
         [0024]    The computer  10  typically includes at least one processing unit  12  (illustrated as “CPU”) coupled to a memory  14  along with several different types of peripheral devices, e.g., a mass storage device  16 , a user interface  18  (including, for example, user input devices and a display), and a network interface  20 . The memory  14  may be comprised of dynamic random access memory (DRAM), static random access memory (SRAM), non-volatile random access memory (NVRAM), persistent memory, flash memory, at least one hard disk drive, and/or another digital storage medium. The mass storage device  16  is typically at least one hard disk drive and may be located externally to the computer  10 , such as in a separate enclosure or in one or more networked computers  21 , one or more networked storage devices  22  (including, for example, a tape drive), and/or one or more other networked devices  23  (including, for example, a server). The computer  10  may communicate with the networked computer, networked storage device  22 , and/or networked device  23  through a network  24 . As illustrated in  FIG. 1 , the computer  10  includes one processing unit  12 , which, in various embodiments, may be a single-thread, multithreaded, multi-core, and/or multi-element processing unit as is well known in the art. In alternative embodiments, the computer  10  may include a plurality of processing units  12  that may include single-thread processing units, multithreaded processing units, multi-core processing units, multi-element processing units, and/or combinations thereof as is well known in the art. Similarly, memory  14  may include one or more levels of data, instruction and/or combination caches, with caches serving an individual processing unit or multiple processing units as is well known in the art. In some embodiments, the computer  10  may also be configured as a member of a distributed computing environment and communicate with other members of that distributed computing environment through the network  24 . 
         [0025]    The memory  14  of the computer  10  may include an operating system  26  to control the primary operation of the computer  10  in a manner that is well known in the art. In a specific embodiment, the operating system  26  may be a Unix-like operating system, such as Linux. The memory  14  may also include at least one application  28 , or other software program, configured to execute in combination with the operating system  26  and perform a task. It will be appreciated by one having ordinary skill in the art that other operating systems may be used, such as Windows, MacOS, or Unix-based operating systems, for example, Red Hat, Debian, Debian GNU/Linux, etc. 
         [0026]    In general, the routines executed to implement the embodiments of the invention, whether implemented as part of an operating system or a specific application, component, algorithm, program, object, module or sequence of instructions, or even a subset thereof, will be referred to herein as “computer program code” or simply “program code.” Program code typically comprises one or more instructions that are resident at various times in memory and storage devices in a computer, and that, when read and executed by at least one processor in a computer, cause that computer to perform the steps necessary to execute steps or elements embodying the various aspects of the invention. Moreover, while the invention has been, and hereinafter will be, described in the context of fully functioning computers and computer systems, those skilled in the art will appreciate that the various embodiments of the invention are capable of being distributed as a program product in a variety of forms, and that the invention applies regardless of the particular type of computer readable media used to actually carry out the invention. Examples of computer readable media include, but are not limited to, recordable type media such as volatile and non-volatile memory devices, floppy and other removable disks, hard disk drives, tape drives, optical disks (e.g., CD-ROM&#39;s, DVD&#39;s, HD-DVD&#39;s, Blu-Ray Discs), among others, and transmission-type media such as digital and analog communications links. 
         [0027]    In addition, various program code described hereinafter may be identified based upon the application or software component within which it is implemented in specific embodiments of the invention. However, it should be appreciated that any particular program nomenclature that follows is merely for convenience; and thus, the invention should not be limited to use solely in any specific application identified and/or implied by such nomenclature. Furthermore, given the typically endless number of manners in which computer programs may be organized into routines, procedures, methods, modules, objects, and the like, as well as the various manners in which program functionality may be allocated among various software layers that are resident within a typical computer (e.g., operating systems, libraries, Application Programming Interfaces [APIs], applications, applets, etc.), it should be appreciated that the invention is not limited to the specific organization and allocation of program functionality described herein. 
         [0028]    Those skilled in the art will recognize that the environment illustrated in  FIG. 1  is not intended to limit the present invention. Indeed, those skilled in the art will recognize that other alternative hardware and/or software environments may be used without departing from the scope of the invention. 
         [0029]    With the computing system described in some detail, an exemplary physical system to be modeled is described. The system to be modeled and simulated includes an etching apparatus described in the plasma model and can be a DF-CCP etching system such as those taught in U.S. Pat. Nos. 5,997,962; 6,426,477; and 6,576,860 each incorporated herein by reference in its entirety; however, the necessary details will be described specifically herein. This system is selected because its nominal complexity renders it suitable for illustration purposes. It should be understood that the principles described herein are applicable to a wide variety of other plasma processing systems. Other examples of plasma processing systems and plasma sources that are described with the aid of computer simulations are described in U.S. Pat. Nos. 7,273,533 and 6,523,493 and U.S. Patent Application Publication No. US2005/0160985, each expressly incorporated by reference herein in its entirety. 
         [0030]      FIG. 2  schematically illustrates the basic components of a representative DF-CCP system  30 , which includes a process chamber  32  formed within a container  34  constructed from conductive material, such as aluminum having an oxidized alumite surface. The container  34  is grounded through a ground line  36 . A support plate  38  formed of an insulating material is arranged on the bottom of the process chamber  32 . A wafer support  40  is arranged on the insulating support plate  38  and is operable to hold the wafer  42  to be etched. The wafer support  40  also acts as a lower electrode. 
         [0031]    The wafer support  40  is made of aluminum having an oxide surface and can include several features, such as heating and chilling elements (not shown). 
         [0032]    The wafer support  40  is surrounded by a baffle plate  44  arranged around the wafer support  40 . An inner portion of the baffle plate  44  is affixed to the wafer support  40 . The baffle plate  44  serves to facilitate the uniform flow of process gas from a processing volume  46  of the process chamber  32  above the wafer  42  to an exhaust volume  48  in the lower portion of the process chamber  32  from which it is exhausted by an exhaust pump  50 . 
         [0033]    A shower head  52  is arranged at the top of the process chamber  32  for introducing a gas from a gas source  54 . The shower head  52  is attached to the ceiling  56  of the process chamber  32  through an insulating support  58  made of alumina and a cooling member  60  made of aluminum. An upper electrode  62  is fixed to the bottom of the shower head  52  and is made of a conductive material, such as single crystalline silicon. The upper electrode  62  serves as the outlet plate of the shower head  52  and includes a plurality of delivery holes  64  through which a gas within the shower head  52  can be delivered into the process chamber  32  above the wafer  42  on the wafer support  40 . The upper electrode  62  should be approximately the same in size as the lower electrode (represented by the wafer support  40 ). 
         [0034]    The process chamber  32  further includes a shield ring  66  near the upper electrode  62  and made of an insulating material such as quartz. A second insulating ring  68  can be made of a fluoride-based resin and positioned between the top of the shield ring  66  and ceiling  56  of the process chamber  32 . Generally, the frequency of the RF for the upper electrode  62  should be much greater than the frequency of the RF to the lower electrode. 
         [0035]    The upper electrode  62  and the lower electrode (represented by the wafer support  40 ) are each respectively connected to individual power sources or RF power supplies  70 ,  72 . Each RF power supply is respectively connected through a matching network  74 ,  76  to the upper electrode  62  and the lower electrode (represented by the wafer support  40 ). Generally, the RF power supply  72  for the lower electrode will operate at a lower frequency, such as 2 MHz (from about 1 MHz to about 5 MHz). The RF power supply  70  for the upper electrode  62  will operate at a higher frequency, such as 100 MHz (from about 10 MHz to about 200 MHz). 
         [0036]    The DF-CCP system  30  includes a control system (not shown) operable to govern the RF power supplies  70 ,  72  and to operate the other sub-systems of the system  30  to perform a plasma etch process on the wafer  42 . 
         [0037]    In operation, the DF-CCP system  30  will generate a capacitively coupled plasma in the processing volume  46  for etching the wafer  42  upon the wafer support  40  according to a first embodiment of the present invention by which, for example, an oxide film such as SiO 2  is etched on the wafer  42 . 
         [0038]    In performing an etching process in the DF-CCP system  30 , a wafer  42  is positioned upon the wafer support  40  and the processing volume  46  is set to a predetermined vacuum. One or more process gases, including Ar, N 2 , CF 4  and others, is supplied from the gas source  54  at a predetermined flow rate. The one or more process gases are exhausted at the exhaust volume  48  by the vacuum exhaust pump  50 . Typical pressures are in the range of from 1 to 20 mTorr (up to about 300 mTorr). 
         [0039]    After the vacuum is established, the upper electrode  62  is supplied with an RF energy having a frequency of 100 MHz and a power of, e.g., 2 kW from the RF power supply  70 . The lower electrode (represented by the wafer support  40 ) is supplied with an RF energy having a frequency of 2 MHz and a power of, e.g., 1 kW from the RF power supply  72 . The RF power applied to the upper electrode  62  generates a plasma between the upper electrode and the wafer support  40 , while the RF power to the lower electrode (represented by the wafer support  40 ) provides a self-bias on the wafer  42  for directing the etchant ions within the plasma onto the wafer  42 . This combination of lower and higher frequency RF energies creates an immediate and stable plasma. 
         [0040]    Within the plasma, electrons are energized by the high frequency RF field. As the high energy electrons diffuse through the plasma, the high energy electrons will undergo non-elastic collisions with molecules of the process gas. These collisions can result in the generation of radicals, cations, and/or anions by the dissociation of the process gas wherein the generating of these species will depend on the particular gas, the energy of the high energy electrons, and the number of collisions occurring between electrons and a particular gas molecule. Additionally, a portion of the dissociated radicals, anions, and cations may undergo recombination into the process gas molecules. 
         [0041]    Cation etchant ions of gas from the plasma are attracted by a negative bias on the wafer support  40  such that the cation etchant ions are incident onto the wafer  42 . As a result, the silicon oxide (SiO 2 ) film on the surface of the wafer  42  can be etched at a uniform etching rate that is determined in part by the composition, density, and distribution of the plasma in the processing volume  46 , and physical parameters such as the geometric, electric, and chemical parameters of the process chamber  32  and the DF-CCP system  30 . 
         [0042]    When the etching is complete, the controller stops the supply of RF energy to the upper electrode  62  and the wafer support  40 , and stops the supply of the process gas into the process chamber  32 , whereupon the wafer  42  can be removed from the process chamber  32 . 
         [0043]    Existing between the plasma and the upper and lower electrode there is a plasma sheath. The plasma sheath properties are usually considered to be a vacuum. 
         [0044]    While the basic components of the DF-CCP system  30  have been described in some detail, it would be understood than additional features or geometries could be included within the simulation model. 
         [0045]    With the DF-CCP system  30  now described with some detail, the simulation according to an embodiment of the present invention for modeling of the distribution and other parameters of the plasma, and the computer method for such modeling, will now be described. The plasma to be simulated is specified in a mathematical coordinate system called the simulation domain. The simulation domain used herein is two-dimensional and is constructed with axial symmetry; however, multiple dimension systems are also possible. The plasma, for example, can be modeled in a manner described by Lee et al., “Modeling electromagnetic effects in capacitive discharges,”  Plasma Sources Sci. Technol.  17, 1-16 (2008), incorporated herein by reference in its entirety. 
         [0046]    The program code to simulate the plasma may be executed as part of, or executed on behalf of, a software suite, application, command, or request. In some embodiments, the program code may be incorporated with, or executed on behalf of, device simulation software. In a specific embodiment, the program code may be incorporated with, or executed on behalf of, a version of COMSOL application/software suite as distributed by The COMSOL Group of Burlington, Mass. In alternative embodiments, the program code may be incorporated with, or executed on behalf of, mathematical software, such as a version of Fluent by ANSYS Corp or Mathematica by Wolfram Research, Inc. 
         [0047]    The transient method of modeling the plasma is shown in  FIG. 3 . This method can, for example, model transient behavior that leads to steady state solutions of various plasma properties. To begin the method, and in accordance with steps  80 ,  82 ,  84  and  86 , configuration settings, initial conditions, boundary conditions, and parameters of the etching system are entered into the start of the simulation program. These initial parameters and settings are utilized within various subsequent steps as shown. 
         [0048]    After the simulation domain and initial parameters are set, the simulation begins at step  88 , wherein the simulation domain is entered. This simulation domain involves solving the transient domain behavior for the time-based variables to arrive at a steady state solution of the plasma. In the DC-CCP system  30 , this involves solving equations that describe the system for each of the time-based variables, which in the case of the system  30  are CCP power that has a first frequency of, for example, 100 mHz and the substrate bias power that has a second frequency of, for example, 2 mHz, in a series of nested loops. 
         [0049]    The transient behavior based on the CCP power at the first frequency is iterated in a first and inner loop  87  over a time interval sufficient for convergence to a steady state solution, which for the 100 mHz frequency is found to be a time interval of from 0 to approximately 1-10 microseconds, in 0.1 nanosecond increments, for a total of about 10,000 to 100,000 iterations. The transient behavior based on the bias power at the 2 MHz frequency is iterated in a second and outer loop  89  over a time interval also deemed sufficient for convergence. Where the first frequency is much larger than the second frequency, this is found to be one complete cycle of 500 nanoseconds, calculated in 10 or 15 degree increments, for example, for twenty iterations of 25 nanosecond. The solution for the system  30  with the two time bases thus requires about 250,000 to 2,500,000 iterations that each involves solving of all of the equations describing the system. 
         [0050]    Execution of the inner loop  89  begins with calculating the RF power deposition within the plasma as in step  90 , using the configuration settings (step  80 ), the initial settings (step  82 ), and the boundary conditions (step  84 ). One manner of calculating the RF power deposition is by an Electro-Magnetic solver for Maxwell Equations in a simple format having 2-D axial symmetry (i.e. electric field with components in the r-z plane and a magnetic field having only an azimuthal component). As such, the RF power deposition can be computed from: 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       ∇ 
                       
                         × 
                         
                           ( 
                           
                             
                               
                                 ( 
                                 
                                   
                                     ɛ 
                                     r 
                                   
                                   - 
                                   
                                     j 
                                      
                                     
                                         
                                     
                                      
                                     
                                       σ 
                                       
                                         ωɛ 
                                         0 
                                       
                                     
                                   
                                 
                                 ) 
                               
                               
                                 - 
                                 1 
                               
                             
                              
                             
                               ∇ 
                               
                                 × 
                                 
                                   H 
                                   ϕ 
                                 
                               
                             
                           
                           ) 
                         
                       
                     
                     - 
                     
                       
                         μ 
                         r 
                       
                        
                       
                         k 
                         0 
                       
                        
                       
                         H 
                         ϕ 
                       
                     
                   
                   = 
                   0 
                 
               
               
                 
                   EQUATION 
                    
                   
                       
                   
                    
                   1 
                 
               
             
           
         
       
     
         [0000]    where ε r  is the relative permittivity in the sheath, j is the complex number index, σ is the conductivity, ω is the RF excitation angular frequency, ε 0  is the permittivity of free space, H φ  is the magnetizing field intensity, μ r  is the permeability, and k 0  is the wave number in vacuum given by k 0 =ω√{square root over (ε 0 μ 0 )}=ω/c 0 , where ω is the angular frequency and c 0  is the speed of light. Solving for H φ  results in the electric field: 
         [0000]      H φ           (E r ,E z )   EQUATION 2 
         [0000]    and the RF power that is deposited into the plasma is calculated as: 
         [0000]    
       
         
           
             
               
                 
                   
                     Q 
                     abs 
                   
                   = 
                   
                     
                       1 
                       2 
                     
                      
                     
                       σ 
                       p 
                     
                      
                     
                       
                          
                         E 
                          
                       
                       2 
                     
                   
                 
               
               
                 
                   EQUATION 
                    
                   
                       
                   
                    
                   3 
                 
               
             
           
         
       
     
         [0000]    where Q abs  is the heat absorbed from the RF field and σ p  is the conductivity of the plasma. 
         [0051]    Once the iteration step on the calculation of the radiofrequency field and power deposition is complete, the gas flow is calculated in step  92 . This can be accomplished from steady state Navier-Stokes equations assuming a low Reynolds number and laminar flow in 2-D axial symmetry geometry. The total stress tensor formulation is given by: 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       ρ 
                        
                       
                         
                           ∂ 
                           
                             u 
                             → 
                           
                         
                         
                           ∂ 
                           t 
                         
                       
                     
                     - 
                     
                       ∇ 
                       
                         · 
                         
                           [ 
                           
                             
                               
                                 - 
                                 p 
                               
                                
                               
                                   
                               
                                
                               
                                 I 
                                 → 
                               
                             
                             + 
                             
                               η 
                               ( 
                               
                                 
                                   ∇ 
                                   
                                     u 
                                     → 
                                   
                                 
                                 + 
                                 
                                   
                                     ( 
                                     
                                       ∇ 
                                       
                                         u 
                                         → 
                                       
                                     
                                     ) 
                                   
                                   T 
                                 
                               
                               ) 
                             
                           
                           ] 
                         
                       
                     
                     + 
                     
                       
                         ρ 
                          
                         
                           ( 
                           
                             
                               u 
                               → 
                             
                             · 
                             ∇ 
                           
                           ) 
                         
                       
                        
                       
                         u 
                         → 
                       
                     
                   
                   = 
                   
                     F 
                     → 
                   
                 
               
               
                 
                   EQUATION 
                    
                   
                       
                   
                    
                   4 
                 
               
             
           
         
       
     
         [0000]    where ρ is the gas density, {right arrow over (u)} is the gas velocity, p is the process chamber pressure, {right arrow over (I)} is the unit diagonal matrix, η is the gas viscosity, and {right arrow over (F)} are the body forces. The mass continuity equation under imcompressible gas assumption (that is with density being constant) will simplify to: 
         [0000]      ∇· {right arrow over (u)}= 0   EQUATION 5 
         [0052]    Turning now to iteration step  94 , the temperature of the gas is calculated in accordance with energy conservation. As noted above, the gas molecules are heated (source term Q) due to the non-elastic collisions between gas molecules and RF field energized electrons. Thus, the gas temperature, T g , is computed through a heat transport by a gas convection and heat conduction module: 
         [0000]      ∇(− k∇T   g )= Q−ρc   p   {right arrow over (u)}∇T   g    EQUATION 6 
         [0000]    where k is Boltzmann constant and c p , ρ, and k are gas heat capacity, density, and thermal conductivity of the gas mixture, respectively. 
         [0053]    With the gas flow and temperature calculated, the local time iterative scheme can begin at step  96 . As noted above, the upper electrode is powered by a 100 MHz RF frequency, which corresponds to a 10 ns time scale. As a result, the iterations to simultaneously solve the equations of step  96  through step  100  occur at 0.1 ns intervals until a transient convergence of the intermediate values of the electron temperature, electron density, and radical densities is achieved. This iteration process can typically continue at about 10 −4  second (or 100 μs) intervals within the internal time of the model. Specifically, the iterative calculations include the following steps and equations. 
         [0054]    In step  96 , the electron temperature (T e ) is computed in accordance with the conservation of energy: 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       
                         ∂ 
                         
                           ∂ 
                           t 
                         
                       
                        
                       
                         ( 
                         
                           
                             3 
                             2 
                           
                            
                           
                             n 
                             e 
                           
                            
                           
                             T 
                             e 
                           
                         
                         ) 
                       
                     
                     + 
                     
                       ∇ 
                       
                         ( 
                         
                           
                             
                               5 
                               2 
                             
                              
                             
                               T 
                               e 
                             
                              
                             
                               Γ 
                               e 
                             
                           
                           - 
                           
                             k 
                              
                             
                                 
                             
                              
                             
                               ∇ 
                               
                                 T 
                                 e 
                               
                             
                           
                         
                         ) 
                       
                     
                   
                   = 
                   
                     
                       Q 
                       abs 
                     
                     - 
                     
                       e 
                        
                       
                           
                       
                        
                       
                         Γ 
                         e 
                       
                        
                       
                         E 
                         a 
                       
                     
                     - 
                     
                       
                         
                           3 
                            
                           m 
                         
                         M 
                       
                        
                       
                         v 
                         m 
                       
                        
                       
                         n 
                         e 
                       
                        
                       
                         T 
                         e 
                       
                     
                     - 
                     
                       
                         ∑ 
                         i 
                       
                        
                       
                         
                           V 
                           i 
                         
                          
                         
                           υ 
                           i 
                         
                          
                         
                             
                         
                          
                         
                           n 
                           e 
                         
                       
                     
                   
                 
               
               
                 
                   EQUATION 
                    
                   
                       
                   
                    
                   7 
                 
               
             
           
         
       
     
         [0000]    where E a  is an ambipolar electric field; Γ e  is electron flux; v i  is the collision frequency of the “i”-th reaction and V i  is its threshold; e is the elementary charge; m is electron mass; and M is total mass of two particles in elastic collision. 
         [0055]    In step  98 , the electron (or ion) density is calculated by an ambipolar-drift diffusion approach, with a transient mass balance equation for “i” type of ions: 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       
                         ∂ 
                         
                           n 
                           i 
                         
                       
                       
                         ∂ 
                         t 
                       
                     
                     + 
                     
                       ∇ 
                       
                         ( 
                         
                           
                             
                               - 
                               
                                 D 
                                 a 
                               
                             
                              
                             
                               ∇ 
                               
                                 n 
                                 i 
                               
                             
                           
                           + 
                           
                             
                               n 
                               i 
                             
                              
                             
                               u 
                               → 
                             
                           
                         
                         ) 
                       
                     
                   
                   = 
                   
                     R 
                     i 
                   
                 
               
               
                 
                   EQUATION 
                    
                   
                       
                   
                    
                   8 
                 
               
             
           
         
       
     
         [0000]    Where n i  is the “i”-th ion and D a  is the ambipolar diffusion coefficient. 
         [0056]    Finally, in step  100 , the neutral species transport within the plasma is calculated from the mass balance equations: 
         [0000]      ∇(− D   k   ∇n   k   +n   k   {right arrow over (u)} )= R   k    EQUATION 9 
         [0000]    where D k  is the diffusion coefficient of the “k”-th particle and n k  is the “k” particle, and where k=all neutral species of the reaction. The generation and recombination of the positive etchant ions is given by: 
         [0000]    
       
         
           
             
               
                 
                   
                     R 
                     k 
                   
                   = 
                   
                     
                       
                         n 
                         e 
                       
                        
                       
                         
                           ∑ 
                           j 
                         
                          
                         
                           
                             n 
                             j 
                           
                            
                           
                             〈 
                             
                               
                                 σ 
                                 jk 
                               
                                
                               
                                 v 
                                 e 
                               
                             
                             〉 
                           
                         
                       
                     
                     + 
                     
                       
                         ∑ 
                         
                           j 
                           , 
                           l 
                         
                       
                        
                       
                         
                           R 
                           jkl 
                         
                          
                         
                           n 
                           j 
                         
                          
                         
                           n 
                           l 
                         
                       
                     
                     - 
                     
                       
                         n 
                         e 
                       
                        
                       
                         n 
                         k 
                       
                        
                       
                         
                           ∑ 
                           j 
                         
                          
                         
                           〈 
                           
                             
                               σ 
                               kj 
                             
                              
                             
                               υ 
                               e 
                             
                           
                           〉 
                         
                       
                     
                     - 
                     
                       
                         n 
                         k 
                       
                        
                       
                         
                           ∑ 
                           
                             j 
                             , 
                             l 
                           
                         
                          
                         
                           
                             R 
                             kjl 
                           
                            
                           
                             n 
                             l 
                           
                         
                       
                     
                   
                 
               
               
                 
                   EQUATION 
                    
                   
                       
                   
                    
                   10 
                 
               
             
           
         
       
     
         [0057]    where n e  is the electron, n j  is the “j” particle, σ jk  and σ kj  are the recombination rates due to collisions with electrons, R jkl  (forward) and R kjl  (backward) are the rate constants due to collision with species n k  or n j , and v e  is the electron velocity. 
         [0058]    As indicated above, iteration steps  90 - 100  are repeated at approximately 0.1 ns intervals until the transient solutions converge. Typically, the iterations involve calculations of steps  96 - 100  for electron temperature, electron density, and radical densities to converge. Convergence usually occurs within an internal model time of between 1 μs and 10 μs, at which time the iteration scheme computation proceeds to loop switch  102 , which involves recalculation of all steps  90 - 100 . 
         [0059]    Throughout step  96  to step  100 , the local constraints and convergence criteria of step  104  can be changed or adjusted to aid in the convergence of the mentioned parameters. 
         [0060]    If convergence does not occur at step  102 , the new iteration sequence will start again at Step  90  through Step  100  under the satisfaction of the global constraints and convergence criteria at step  105 . The process begins again at step  90  with these new criteria. 
         [0061]    When the electron density for the first frequency has converged, then the RF cycle at the second frequency is computed in the outer loop  89  at phase steps that are less than 5% of the RF cycle duration, that is phase modulated by 10° to 15° steps, with the process of steps  88  through  106  are repeated for each iteration of loop  89 . Typically, steps  90  through  102  can take about 4 hours to about 6 hours to complete the coupled iteration scheme computation using a dual core PC at 2-3 GHz. Thus, a full transient solution for a dual frequency system can take nearly 100+ hours. 
         [0062]    To decrease the computational time needed for achieving a solution, the transient model is converted to a stochastic model as generally shown in  FIG. 4A  and in accordance with the present invention. 
         [0063]    Turning now to Step  108  of  FIG. 4A , the method begins with defining the simulation domain according to one embodiment of the present invention. The simulation domain defining the characteristics of the etching system will include time dependent variables that are dependent upon the chemicals, materials, or geometries of the plasma. Often, the simulation domain may contain at least one or more time dependent variables, or sub-domains for differing chemicals, materials, and/or geometries. That is, for each sub-domain, the time dependence of the variable in each point within the sub-domain should be a homogenous transient formulation, i.e. have similar time dependence differing only in phase for various points. Otherwise, the simulation domain is a heterogenous transient solution and should be divided into two or more sub-domains of homogeneous transient formulations. 
         [0064]    For example, as illustrated in  FIG. 5 , the simulation domain D 2    110  can include points (x 1 ,z 1 ) and (x 2 ,z 2 ) shown to differ in their time dependence of a particular variable. Thus, points (x 1 ,z 1 ) and (x 2 ,z 2 ) will exist in first and second sub-domains (D 2   1 , D 2   2 )  112 ,  114 , respectively. While only two sub-domains D 2   1  and D 2   2    112 ,  114  within the simulation domain D 2    110  are shown, it would be understood that more sub-domains are possible. 
         [0065]    A time dependent function, Ψ, for the time dependent variable of each sub-domain can be determined by empirical data or by analytical derivation and are represented by the graphs Ψ versus time within  FIG. 5 . 
         [0066]    Returning again to  FIG. 4A , and specifically to Step  116  where a probability function Θ for each time dependent function, Ψ, is determined, wherein the probability function Θ is a statistical distribution of the values of Ψ. The probability function can be determined from the known analytical expression of the transient dependence of the particular variable by measured or empirically established transient dependence, hypothetically derived, or measured directly. As seen in  FIG. 5 , each Ψ of a separate sub-domain may have a different statistical distribution shown in graphs Θ versus Ψ. 
         [0067]    Continuing with  FIG. 4A , and step  118 , the probability function Θ is normalized to ensure that an occurrence of all transient values corresponds to the probability equal to 1. In other words, the probability function is multiplied by a normalizing constant, κ according to: 
         [0000]      ∫ Ψ     min     Ψ     max   κΘ(Ψ) dΨ= 1   EQUATION 11 
         [0068]    This will ensure that the area under the Θ versus Ψ curve (see  FIG. 6 ) should be equal to 1. The normalized probability function will then be given by: 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       
                         Θ 
                         norm 
                       
                        
                       
                         ( 
                         Ψ 
                         ) 
                       
                     
                     = 
                     
                       κΘ 
                        
                       
                         ( 
                         Ψ 
                         ) 
                       
                     
                   
                   , 
                   
                     
                       where 
                        
                       
                           
                       
                        
                       κ 
                     
                     = 
                     
                       1 
                       
                         
                           ∫ 
                           
                             Ψ 
                             min 
                           
                           
                             Ψ 
                             max 
                           
                         
                          
                         
                           
                             Θ 
                              
                             
                               ( 
                               Ψ 
                               ) 
                             
                           
                            
                           
                              
                             Ψ 
                           
                         
                       
                     
                   
                 
               
               
                 
                   EQUATION 
                    
                   
                       
                   
                    
                   12 
                 
               
             
           
         
       
     
         [0069]    With the probability function now normalized, the method proceeds to step  120  of  FIG. 4A  with reference made to  FIG. 6 . In step  120 , the interval of the function Ψ values, that is &lt;Ψ min ,Ψ max &gt; is subdivided into “n” subintervals &lt;Ψ min :Ψ max &gt; nks , wherein the size of each subinterval will be given by: 
         [0000]    
       
         
           
             
               
                 
                   
                     ΔΨ 
                      
                     
                       ( 
                       
                         x 
                         , 
                         z 
                       
                       ) 
                     
                   
                   = 
                   
                     
                       
                         
                           Ψ 
                           max 
                         
                          
                         
                           ( 
                           
                             x 
                             , 
                             z 
                           
                           ) 
                         
                       
                       - 
                       
                         
                           Ψ 
                           min 
                         
                          
                         
                           ( 
                           
                             x 
                             , 
                             z 
                           
                           ) 
                         
                       
                     
                     
                       n 
                       ks 
                     
                   
                 
               
               
                 
                   EQUATION 
                    
                   
                       
                   
                    
                   13 
                 
               
             
           
         
       
     
         [0000]    And the characteristic value,  Ψ   b , is determined for each n ks  sub-interval. The  Ψ   b  can be mean value of Ψ values between &lt;Ψ min :Ψ max &gt;; however, other values may also be appropriate. Additionally, the subintervals should be sufficiently small to minimize sampling error. This is illustrated in  FIG. 6 , wherein the characteristic value  Ψ  for subintervals  Ψ   b−1 ,  Ψ   b , and  Ψ   b+1  are shown. 
         [0070]    Turning again to  FIG. 4A , and specifically to Step  122 , wherein the characteristic values  Ψ   b  are arranged into corresponding strings, B n , wherein the length of Bn is given by: 
         [0000]      Δ b=INT[n   NS   n   MS Θ norm [  Ψ   b ( x,z )]]  EQUATION 14 
         [0000]    where n NS  and n MS  are the number of points within the simulation domain  110  in the x- and z-axes, respectively, and each member within the string will equal the value of  Ψ   b . 
         [0071]    In step  124 , all strings of B n  are arranged into an array Ω corresponding to each sub-domain, as shown in  FIG. 7 . 
         [0072]    In step  126  of  FIG. 4A  and further in  FIG. 7 , the values within each Ω are reassigned randomly into the points of the corresponding sub-domain. Once all of the points of the original simulation domain are randomly reassigned values from the respective Ω, then the resultant array is a stochastic simulation array  128  and is ready to undergo simulation, as described in  FIG. 4B . 
         [0073]      FIG. 4B  illustrates the simulation of the DF-CCP device with the stochastic simulation array  128  created from step  130  ( FIGS. 4A and 4B ). While the majority of the simulation method is similar to what was detailed previously with  FIG. 3 , distinct differences will be noted. Most notably is the use of the stochastic simulation array eliminates the need for the super-loop or outer loop  89  with a phase-modulation at the second frequency in the iterative scheme. 
         [0074]    Described herein is a method and an algorithm for converting a transient to spatial stochastic formulation for enhancing the computation of plasma fluid models. The computer method may be encoded on a computer usable medium and provided in the form of a computer program product. The system of  FIG. 1 , when programmed with the computer program product, is a simulation apparatus that includes a calculating unit programmed to calculate a probability function as described in connection with  FIG. 4A  and a simulation unit programmed to perform a simulation of transient plasma sheath behavior as described in connection with  FIG. 4B . 
         [0075]    While the invention has been illustrated by a description of various embodiments and while these embodiments have been described in considerable detail, it is not the intention of the applicant to restrict or in any way limit the scope of the appended claims to such detail. Additional advantages and modifications will readily appear to those skilled in the art. Thus, the invention in its broader aspects is therefore not limited to the specific details, representative apparatus and method, and illustrative example shown and described. Accordingly, departures may be made from such details without departing from the spirit or scope of applicants&#39; general inventive concept.