Abstract:
Characteristics of a plasma contained in a reaction chamber of a plasma reactor are determined by first computing plasma characteristics for each of a plurality of cross-sections of the reaction chamber, and then generating a generalized model of the plasma from the computed plasma characteristics for the plurality of cross-sections, for example, by averaging the computed plasma characteristics for the cross-sections. The plasma reactor may comprise a plurality of magnets that move with respect to the reaction chamber, such as in a dipole ring magnet (DRM) plasma reactor, and each of the plurality of cross-sections may include an axis of rotation about which the magnets rotate. Plasma characteristics for each the cross-sections of the reaction chamber may be computed by computing electron density and temperature using a Monte Carlo computational procedure and computing ion and neutral species transmission phenomena from a plasma dynamics simulation, e.g., by computing solutions to a continuity equation and Poisson&#39;s equation for the ion and neutral species. A static magnetic field generated by the moving magnets may be determined, and the plasma characteristics for each of the plurality of cross-sections may be from the determined static magnetic field, shape information for the reaction chamber, and plasma collision reaction data. The generalized model may be used, for example, to estimate an etching rate for a wafer positioned in the chamber.

Description:
RELATED APPLICATION  
       [0001]    This application is related to Korean Application No. 2001-167, filed Jan. 3, 2001, the disclosure of which is hereby incorporated herein by reference. 
     
    
     
       BACKGROUND OF THE INVENTION  
         [0002]    The present invention relates to methods, apparatus and computer program products for simulating plasma behavior in a plasma reactor apparatus, such as those widely used for manufacturing semiconductor devices.  
           [0003]    In 1996, World Semiconductor Statistics (WSTS) showed that plasma-related equipment accounted for 40% of all semiconductor manufacturing equipment sales. Plasma processes are extensively used for deposition, ion implantation, cleaning, and etching. The use of plasma processes in manufacturing semiconductor devices is expected to increase.  
           [0004]    Plasma etching processes used in manufacturing highly integrated semiconductor devices generally require precise control to meet requirements such as uniformity, selectivity ratio and anisotropy. Thus, setting up a mass production process using plasma etching techniques can be costly and time-consuming.  
           [0005]    Such cost and time may be reduced by simulating plasma behavior. In particular, process development generally requires understanding of surface reaction and other phenomena associated with the plasma processing. Thus, plasma modeling and simulation can be valuable.  
           [0006]    [0006]FIG. 1 is a flowchart of a conventional simulation method for inductively coupled plasma (ICP) equipment. Referring to FIG. 1, plasma reactor shape and process conditions (block  2 ) and data on plasma collision reaction (block  4 ) are provided. A plasma simulation (block  6 ) comprises three operations: a module that determines the electromagnetic field (block  8 ), a module that calculates electron density and temperature using a Monte Carlo technique (block  10 ), and a module that determines transmission phenomena of chemical species (block  12 ). These three operations are repeated until they converge to a result. This simulation results in estimates for plasma characteristics (block  14 ), such as electromagnetic field distribution, electron density and temperature, ion and neutral species distribution directly involved in surface reaction, and flux incident onto a wafer surface in a plasma reactor, all of which can affect etching processes. However, such simulations typically employ a three-dimensional calculation that can take several days or longer. Therefore, it may be impractical to apply such a simulation approach in the development of a real plasma process.  
           [0007]    Plasma etching processes used in manufacturing semiconductor devices typically use dipole ring magnet (DRM) plasma equipment. Typical DRM plasma equipment implements a magnetically enhanced reactive ion etching (MERIE) method using a complex structure that includes several (e.g.,  20 ) permanent magnets having different magnetic forces and fluxes that rotate around a plasma reaction chamber at speeds on the order of 20 revolutions per minute (rpm) (See “A New High-Density Plasma Etching System Using a Dipole-Ring Magnet”, JJAP, pp. 6274-6278, 1995).  
           [0008]    A plasma having external magnetic fields applied thereto may be simulated using a conventional 3-dimensional calculation method (See “A three-dimensional model for inductively coupled plasma etching reactors: Azimuthal symmetry, coil properties, and comparison to experiments”, JAP, pp. 1337-1344, 1996). However, as discussed above, a conventional 3-dimensional simulation method may require a calculation time of several days or more. Therefore, it may be impractical to apply such a conventional 3-dimensional simulation method to the development of a real process using a structure such as that found in DRM plasma equipment.  
         SUMMARY OF THE INVENTION  
         [0009]    According to embodiments of the present invention, characteristics of a plasma contained in a reaction chamber of a plasma reactor are determined. Plasma characteristics for each of a plurality of cross-sections of the reaction chamber are first determined, and then a generalized model of the plasma is generated from the computed plasma characteristics for the plurality of cross-sections. For example, the plasma reactor may comprise a plurality of magnets that move with respect to the reaction chamber, such as in a dipole ring magnet (DRM) plasma reactor, and each of the plurality of cross-sections may include an axis of rotation about which the magnets rotate.  
           [0010]    In some embodiments of the present invention, computing plasma characteristics for each of a plurality of cross-sections of the reaction chamber comprises computing electron density and temperature for a cross-section using an iterative Monte Carlo computational procedure and computing ion and neutral species transmission phenomena for the cross-section from a plasma dynamics simulation. Computing ion and neutral species transmission phenomena for the cross-section from a plasma dynamics simulation may comprise computing solutions to a continuity equation and Poisson&#39;s equation for the ion and neutral species. Prior to these computations, a static magnetic field generated by the moving magnets may be determined, and the computation of plasma characteristics for each of the plurality of cross-sections of the reaction chamber may comprise computing the plasma characteristics for each of the plurality of cross-sections from the determined static magnetic field, shape information for the reaction chamber, and plasma collision reaction data. Generating a generalized model of the plasma from the computed plasma characteristics for the plurality of cross-sections may comprise computing at least one of an electron density distribution, a temperature distribution, a distribution of ion species, a distribution of neutral species, and a flux incidence, e.g., by averaging the results of the computations performed for the two-dimensional cross-sections. The generalized model may be used, for example, to estimate an etching rate for a wafer positioned in the chamber.  
           [0011]    The present invention may be embodied as methods, apparatus and computer program products. 
       
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS  
       [0012]    [0012]FIG. 1 is a flowchart of a conventional simulation method for inductively coupled plasma (ICP) equipment.  
         [0013]    [0013]FIGS. 2 and 3 gare drawings illustrating a dipole ring magnet (DRM) plasma reactor apparatus.  
         [0014]    [0014]FIGS. 4A and 4B are flowcharts illustrating apparatus and operations for simulating plasma behavior according to embodiments of the present invention.  
         [0015]    [0015]FIG. 5 illustrates a magnetic field induced by magnets of a plasma reactor apparatus.  
         [0016]    [0016]FIG. 6 is a graph illustrating simulated etch rate distributions for a silicon oxide layer obtained from a simulation according to embodiments of the present invention.  
         [0017]    [0017]FIG. 7 is a graph comparing measured plasma density and simulated plasma density as generated by a plasma simulation according to embodiments of the present invention.  
         [0018]    [0018]FIG. 8A is a graph illustrating an etch rate distribution for a silicon oxide layer as a function of etch gas composition ratio estimated according to embodiments of the present invention.  
         [0019]    [0019]FIG. 8B is a graph illustrating an etch rate distribution for a silicon nitride layer as a function of etch gas composition ratio estimated according to embodiments of the present invention. 
     
    
     DETAILED DESCRIPTION  
       [0020]    [0020]FIGS. 2 and 3 gare, respectively, a plane view showing the arrangement of permanent magnets in a DRM plasma apparatus and a cross-sectional view of a DRM plasma apparatus. Referring to FIGS. 2 and 3, a DRM plasma apparatus  100  implements a magnetically enhanced reactive ion etching (MERIE) method and has a structure including about 20 permanent magnets  102  having different magnetic forces and fluxes that rotate around a plasma reaction chamber  101 . The permanent magnets  102  rotate around an axis of rotation  106 . As shown in FIG. 2, magnetic fields  103  of the permanent magnets  102  are arranged in different directions, and form a composite magnetic field  111  in the reaction chamber  101 . The permanent magnets  102  may be differently arranged depending on the type of equipment used. For example, in FIG. 2, the permanent magnets  102  are regularly spaced, while in FIG. 5, the permanent magnets  102  are irregularly spaced.  
         [0021]    The permanent magnets  102  induce a magnetic field  111  that is approximately static, i.e., that is minimally affected by the state of a plasma  104  in the plasma reaction chamber  101 . The time required for stabilizing the plasma  104  in the plasma reaction chamber  101  typically is on the order of hundreds of microseconds or less. A wafer W positioned in the plasma reaction chamber  100  is supported by a chuck C. An electrode  108  is connected to radio frequency power source  110 .  
         [0022]    [0022]FIGS. 4A and 4B are flowchart illustrations of methods, apparatus (systems) and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations, and combinations of blocks, can be implemented by computer program instructions. These computer program instructions may be provided to a processor of a general purpose computer, special purpose computer, or other programmable data processing apparatus, such as mainframe computer, high-performance computer workstation, or parallel-processing system, to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create structures for implementing the functions specified in the block diagram and/or flowchart block or blocks. These computer program instructions may also be stored in a computer-readable memory that can direct a computer or other programmable data processing apparatus to function in a particular manner, such that the instructions stored in the computer-readable memory produce an article of manufacture including instructions which implement the function specified in the block diagram and/or flowchart block or blocks. The computer program instructions may also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer implemented process or method such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the block diagram and/or flowchart block or blocks. Accordingly, the flowcharts of FIGS. 4A and 4B support methods, apparatus and computer program products for performing operations described therein.  
         [0023]    In greater detail, FIGS. 4A and 4B illustrate operations and apparatus for simulating behavior of a plasma in a plasma reactor apparatus such as that illustrated in FIGS. 2 and 3. The configuration of the plasma reactor  100  and process conditions are input (block  20 ). In particular, shape information for the plasma reactor  100 , such as the size of the plasma reaction chamber  101  and the position of the magnets  102 , is input into a simulation program, along with process conditions such as power, pressure, and gas composition ratio. Plasma collision reaction data are also input into the simulation program (block  22 ). Plasma collision reaction data may include a reaction rate constant of the following collision reaction equation, and may be represented by a function such as electron temperature:  
         e — Ar         Ar + +e — +e —   
         e — Cl 2           Cl+Cl+e —   
         [0024]    A 3-dimensional magnetic field induced by the permanent magnet  102  is computed using, for example, commercially-available software (block  24 ). For example, a commercial finite element analysis tool such as Vector Fields may be used to determine the 3-dimensional magnetic field. The magnetic field induced by the permanent magnet  102  is an approximately static magnetic field, which typically is minimally affected by the state of the plasma  104  gin the plasma reactor  100 . Therefore, it is possible to calculate the magnetic field induced by the permanent magnet  102  apart from effects of the plasma  104 .  
         [0025]    Electron density and temperature are computed by a Monte Carlo method (block  30 ) and transmission phenomena of ion and neutral species (block  32 ) are determined until convergence is achieved (block  34 ). In FIG. 4A, electron density and temperature are first calculated by the Monte Carlo method and then the transmission phenomena of ion and neutral species are interpreted. As shown in FIG. 4B, determination of transmission phenomena of ion and neutral species (block  30 ′) may occur before calculation of electron density and temperature (block  32 ′).  
         [0026]    The determination of electron density and temperature (blocks  32 ,  32 ′) and the determination of transmission phenomena of ion and neutral species (blocks  30 ,  30 ′) are performed for 2-dimensional cross-sections of the reaction chamber  101  in a characteristic magnetic field direction. In detail, convergence values ( 36 ) are obtained for each of a plurality of 2-dimensional cross-sections including the axis  106  of rotation.  
         [0027]    The time required for stabilizing the plasma  104  for the given 2-dimensional static magnetic field distribution symmetrical to the axis is on the order of hundreds of microseconds (μs) or less. The permanent magnets  102  typically rotate around the plasma reaction chamber  101  at a speed of about 20rpm, which leads the magnetic field distribution to change depending on time. It is possible to 2-dimensionally sample the cross-sectional magnetic field distribution including the axis  106  in the characteristic magnetic filed direction. The time required for the simulation typically depends on the nature of the plasma  104  and the process conditions. For example, in the case of argon (Ar) plasma under typical process conditions, the simulation may take about one hour.  
         [0028]    According to embodiments of the present invention, calculation of electron density and temperature by Monte Carlo simulation and determination of transmission phenomena of ion and neutral species are performed at a plurality of 2-dimensional cross-sections including an axis for cross-sectional magnetic field distribution in a characteristic magnetic field direction. The convergence values for the sections may be averaged to generate a generalized model of the plasma  104 , for example, electron density and temperature in the plasma reaction chamber  101 , the distribution of ion and neutral species involved in surface reaction at the wafer W, and flux incident onto a major surface of the wafer W.  
         [0029]    Examples of methods for calculating electron density and temperature and of determining transmission phenomena of ion and neutral species will now be described. Collision probabilities of electron-ion, electron-neutral molecule/atom may be calculated and kept as a probability array. The collision frequency v ij  may be expressed as:  
               v   ij     =         (       2        ɛ   i         m   e       )       1   /   2            σ   ij          N   j               (   1   )                               
 
         [0030]    where σ ij  presents electron impact cross-section in I-energy and j-process and N j  presents the density of collision partner in j-process. Ε I  and M e  denote energy and electron mass, respectively. As a result, a probability array P ij  may be expressed by the formula:  
               P   ij     =       [         ∑     l   =   1     j          v   ij       +     (       v   m     -     v   i       )       ]     /     v   m               (   2   )                               
 
         [0031]    where v i  and v m  are respectively expressed by formulas (3) and (4):  
               v   i     =       ∑     l   =   1       l                 max            v   il               (   3   )                               
 
         v m =max(v i )  (4) 
         [0032]    where I max  represents the total number of processes.  
         [0033]    The initial rate and position of an electron, respectively, may be extracted from a Maxwell distribution and a random distribution, and the trajectory of respective pseudo electrons may be separately tracked. A time step Δtl for determination of particle motion may be expressed by formula (5):  
         Δt l =min(00.1τ rf , 00.1τ ECR , t cl −t l )  (5) 
         [0034]    where τ l  is the time until the trajectory of particle l is updated, τ rf  is the radio frequency period, τ ECR  is the local electron cyclotron period, and t ct  is the time until next collision, namely, t ct =t to +V m   −1 1n(R), RΕ[0,1] (where t lO  represents the initial time). A specific process is also selected among several possible processes using a random number generator.  
         [0035]    This Monte Carlo iteration continues for about 20-50 RF cycles (about 3 μs), and an electron impact source function is obtained from a time-averaged electron energy distribution function f (Ε, r, and z) and formula (6):  
               S   ij     =         n   e          (     r   ,   z     )            δ   ij            N   ij          (     r   ,   z     )       ×       ∫   0   ∞            f        (     ɛ   ,   r   ,   z     )            (       2   ɛ       m   e       )            σ   ij          (   ɛ   )                          ɛ                   (   6   )                               
 
         [0036]    where n e  and N ij  respectively represent electron density and collision partner density in I-energy and j-process calculated from a transmission phenomenon interpretation module just before repeating, and τ represents energy. Also, if process ij is a source of j species, δ ij  is +1, gand if process ij is a loss, δ ij  is −1.  
         [0037]    The transmission phenomena determination may involve solving a continuity equation and Poisson&#39;s equation for all ion and neutral species, as expressed by formulas (7) and (8):  
                 δ                   N   j         δ                 t       =       ∇     (         μ   j          q   j          N   j            E   →     s       -       D   j          ∇     N   j           )       +       (       δ                   N   j         δ                 t       )     c               (   7   )                 ∇     ·       E   →     s         =       -       ∇   2        Φ       =     ρ     ɛ   0                 (   8   )                               
 
         [0038]    where μ j , D j , q j , p, (δN j /δt) c , E s , Φ, and Ε O , are the mobility of j-species, diffusion coefficient of j-species, charge of j-species, charge density, density variation by all collisions, electric field, electrostatic potential, and dielectric constant of a vacuum state, respectively. (δN j /δt) c  includes contribution by heavy particles as well as contribution from S j  (generation rate in coordinates r and z), which are not distinguished in Formulas (7) and (8).  
         [0039]    The continuity equation may exhibit a problem where the Knudsen number λ/L (λ is an averaged free path and L is the length of a reactor) is increased to greater than 0.1 with less than 100 mTorr pressure, diffusion velocity may get faster than the thermal velocity (V th ) of respective species during drift-diffusion. Consequently, to prevent this phenomenon, diffusion coefficient and particle mobility may limited, as expressed by formulas (9) and (10):  
         D j =min(V th L, D j )  (9)                μ   j     =       eD   j       kT   j               (   10   )                                 
         [0040]    where e, k, and T j  are the charge of the electron, Boltzman&#39;s constant, and the temperature of j species, respectively.  
         [0041]    A general plasma dynamics simulation may separately solve Poisson&#39;s equation and the continuity equation. However, when the general plasma dynamics simulation solves these equations simultaneously, it may exhibit a time-step problem. In a case where a transport equation is obtained from explicit differencing, time-step may be limited by a courant limit, as expressed by formula (11):  
               Δ                   t   c       ≤     min        (         Δ                 r         μ   j          E   r         ,       Δ                 z         μ   j          E   z           )               (   11   )                               
 
         [0042]    where Δr and Δz represent spacial mesh sizes, E r  represents an electromagnetic field in a r direction, and E z  represents an electromagnetic field in a z direction. In the case of obtaining an implicit solution, the time-step may be theoretically much larger than the courant limit. However, Poisson&#39;s equation is typically solved by an explicit method regardless of the transport equation. This is why charge density in the current step may be required for updating a potential in a subsequent step, as shown in formula (12):  
                 ∇   2          Φ        (     t   +     Δ                 t       )         =     -       ρ        (   t   )         ɛ   0                 (   12   )                               
 
         [0043]    In this case, the maximum time-step is shorter than a dielectric relaxation time so that the electromagnetic field changes the sign of the time-step, as expressed by formula (13):  
               Δ                   t   d       =       ɛ   0     σ             (   13   )                               
 
         [0044]    where σ represents plasma conductivity. Estimating a dielectric relaxation time value from some calculations, in the case of plasma with low pressure and high density, σ is about 0.1˜1 (Ωcm) −1 . Therefore, Δt d , which is about 10 −13 ˜10 −12  seconds, is much shorter than the courant limit.  
         [0045]    To solve such a short-time-step problem, a semi-implicit differencing type technique may be used in determining transmission phenomena, as expressed by formula (14):  
                 ∇   2          Φ        (     t   +     Δ                 t       )         =     -       1     ɛ   0            [       ρ        (   t   )       +     Δ                 t               ρ        (   t   )              t           ]                 (   14   )                               
 
         [0046]    where a time-derivative of charge density includes only a transport term. The final equation to be solved for which the transport term of formula (7) is expressed by formula (15):  
                   ∇   2          Φ        (     t   +     Δ                 t       )         +       1     ɛ   0          Δ                 t          ∑   i          e                   q   i          μ   i     ×     [         ∇     N   i            ∇     Φ        (     t   +     Δ                 t       )           +       N   i            ∇   2          Φ        (     t   +     Δ                 t       )             ]             =       -       ρ        (   t   )         ɛ   0         -       1     ɛ   0          Δ                 t                 e          ∑   i            q   i          (       ∇     D   i            ∇     N   i         )                     (   15   )                               
 
         [0047]    The formula (15) may be solved by a succession of relaxation (SOR) method, where the optimized SOR parameter is 1.8≦α≦1.9. According to the above semi-implicit technique, the time-step may be about 100˜1000 times lager than Δtd or as large as the courant limit.  
         [0048]    Some have reported that, in the case where the time-step gets larger by the above method, the difference in accuracy is mostly within a few percentage points. However, it has been confirmed that plasma potential and plasma density are about 30% different from an absolute value.  
         [0049]    In solving Poisson&#39;s equation, boundary conditions may depend on whether the surface of a reactor or a substrate is metal or dielectric. In the case where the surface is metal, the surface is grounded or is determined by an external potential. In the case where the surface is a dielectric, potential 40 of the portion in contact with plasma is expressed by formula (16):  
               Φ   0     =       [       Φ   l     +     Δ                   z        [         σ   s     /     ɛ   0       +       ɛ   d            Φ   b     /     (       ɛ   0        L     )           ]               1   +     Δ                 z                     ɛ   d     /     (       ɛ   0        L     )                     (   16   )                               
 
         [0050]    where Φ l  represents the plasma potential at the first mesh on the surface, E d  represents permissivity of the dielectric, L represents the thickness of dielectric, and Φ b  represents plasma potential opposite to the surface. σ s  represents surface charge density, as obtained from formula (17):  
               σ   s     =       ∑   j          ∫     e                   q   j          ∅   j             t                   (   17   )                               
 
         [0051]    where  j =(q j , u j N j {right arrow over (E)}−D j ∇N j )·{circumflex over (n)} represents flux which reaches the surface.  
         [0052]    Convergence velocity generally depends on how close the initial guess for species density approaches an actual value. The time for convergence is typically about 10˜100 μs. Thus, in the case of a low initial guess value, a substantial amount of calculation time may be required. For example, in a case where the time-step is 1 ns (10 −9  s), the time-step may require 10 5  cycles to approach up to 100 μs. Also, 0.025 seconds are typically required per cycle when processed using a Silicon Graphics® Onyx® workstation. Thus, it may take 7 hours to determine just the transmission phemonema.  
         [0053]    To reduce such computation time, an acceleration technique that improves the initial guess using prior results before determination of transmission phenomena may be used. Such an acceleration technique scales up or scales down dN/dt calculated by the initial time-step of about 1·100 ns to about 1000·2000 times to increase the effective time-step to 1000·2000 times, as expressed by formula (18):  
                 N   j          (     r   ,   z   ,     t   +     Δ                 t         )       =         N   j          (     r   ,   z   ,   t     )       +       γ        (            N   j            t       )          Δ                 t               (   18   )                               
 
         [0054]    That is, the parameter Υ for determining acceleration is increased to about 1000·2000, which increases the effective time-step. In this case, converged results may be obtained in about 100·1000 cycles (the maximum number of cycles is generally about 500 in the input step). The density of negatively-charged species may be re-normalized to be equal to the sum of negatively-charge species and positively-charged species, which can solve the charge neutrality problem.  
         [0055]    [0055]FIG. 5 shows an electromagnetic field distribution on a wafer W induced by a permanent magnet. The magnetic field is formed on the wafer W in the plasma reactor  100  by the permanent magnets  102  which are asymmetrically arranged and rotate around the plasma reaction chamber  101 . The magnetic field is a substantially static magnetic field, magnetic flux density of which varies with location on the wafer. For example, on the wafer, magnetic flux density at the point A is about 180 Gauss, the magnetic flux density at the point B is about 120 Gauss, and the magnetic flux density at the point C is about 60 Gauss. According to embodiments of the invention described above, a 2-dimensional plasma simulation is performed for cross-sectional magnetic field distribution in characteristic magnetic field directions, for example, directions I, II, and III. From these results, a generalized plasma behavior model can be generated.  
         [0056]    [0056]FIG. 6 is a graph illustrating an etch rate distribution of a silicon oxide (SiO 2 ) layer calculated according to embodiments of the present invention, for cross-sectional magnetic field distributions comprising 3 sections and 8 sections, respectively. The 2-dimensional plasma simulation was performed assuming a pressure of 25 mTorr, a RF power of 1200 W, CHF 3  flux of 150 sccm, CO flux of 50 sccm, and O 2  flux of 10 sccm. The location on the wafer denotes the distance from the center of the wafer.  
         [0057]    As shown in FIG. 6, the calculated etch rate of the silicon oxide layer at the wafer center is about 1600 Å/min for both the simulation using 3 cross-sections and the simulation using 8 cross-sections. Etch rates of the silicon oxide layer at about 6 cm point from the wafer center are about 1800 Å/min for the simulation at 3 sections and about 1750 Å/min for the simulation value at 8 sections.  
         [0058]    [0058]FIG. 7 is a graph illustrating actual measured plasma density as a function of power and plasma density obtained from a plasma simulation according to embodiments of the present invention using 3 two-dimensional cross-sections. Here, argon (Ar) plasma is used, Ar flux is 200 sccm and pressure is 40 mTorr. As can be seen in FIG. 7, there is close agreement of the measured values and the calculated values from the simulation.  
         [0059]    Table 1 shows simulated and measured values for etch rates of a silicon oxide (SiO 2 ) layer and a silicon nitride (Si 3 N 4 ) layer for various in process conditions. Simulations and experiments were performed for the SiO 2  layer and the Si 3 N 4  layers using a varying etch gas composition ratio and power at a pressure of 35 mTorr.  
                                                                                             TABLE 1                               Etch Rate (Å/min)                SiO 2     Si 3 N 4              Process Conditions   Calculation   Experimental       Calculation   Experimental           (Power/CHF 3 (sccm)/   Value   Value   Errors   Value   Value   Errors       CO(sccm)/O 2 (sccm)   (Å/min)   (Å/min)   (%)   (Å/min)   (Å/min)   (%)                    1500 W/31/150/10   217   2019   4.87   1850   1838   0.78       1500 W/35/150/6   2454   2560   −4.18   2133   2036   4.72       1500 W/39/150/2   2755   2737   0.65   2362   2262   4.39       1200 W/35/150/6   2156   2294   −6.02   1846   1924   −4.04       1800 W/35/15O/6   2665   2726   −2.23   2314   2191   5.60                  
 
         [0060]    As demonstrated in Table 1, the etch rate increases with the fraction of CHF 3 , which is due to the increase in the entire flux and radical flux. An increase in power appears to cause the etch rate to increase. As shown in Table 1, the simulated etch rates of both silicon oxide and silicon nitride layers show good agreement with experimental data, with less than a 6% error.  
         [0061]    [0061]FIG. 8A and 8B are graphs illustrating etch rate distribution of silicon oxide and silicon nitride layers for varying etch gas composition ratio. FIGS. 8A and 8B denote etch rate distributions of the silicon oxide and silicon nitride layers, respectively, with ‘Sim’ denoting calculated values obtained from simulation and ‘Exp’ denoting experimental values obtained by measurement of an actual process. Gas composition ratio terms 31/150/10 are CHF 3  flux (sccm)/CO flux (sccm)/O 2  flux (sccm), respectively. Position on the wafer means the distance from the wafer center. In the case of the silicon oxide layer, the calculated value obtained from simulation at the wafer edge is lower than the experimental value obtained from the actual process. However, simulated etch rates show less than a 6% error. In the case of the silicon nitride layer, the etch rate increases toward the wafer edge, has a maximum value at a predetermined position, and decreases at the edge. Simulated etch rates also show less than a 6% error.  
         [0062]    According to embodiments of the present invention, plasma behavior is simulated using calculations for 2-dimensional cross-sections including an axis of magnet rotation in a characteristic magnetic field direction. As a result, the time needed for simulation can be substantially reduced and the plasma characteristics can be precisely estimated. For example, plasma simulation for a DRM plasma reactor may be performed in a relatively short time, for example, within about 1˜2 hours. The plasma characteristics, such as plasma density and temperature, density distributions of respective chemical species, and flux distribution incident onto the wafer, can be precisely estimated. Based on the plasma characteristics, etch and deposition rates can be estimated. The method, apparatus and computer program products of the present invention can be effectively used for process development and process optimization.