Assessment of conservative weighting scheme in simulating
chemical vapour deposition with trace species
J.-S. Wu
∗;†;‡, W.-J. Hsiao
§, Y.-Y. Lian
§and K.-C. Tseng
§Department of Mechanical Engineering; National Chiao-Tung University; Hsinchu 30050; Taiwan
SUMMARY
Low-pressure or ultra-high vacuum chemical vapour deposition often involves important trace species in both gas-phase and surface reactions. The conservative weighting scheme (J. Thermophys. Heat Transfer 1996; 10(4):579) has been used to deal with the trace species often involved in some non-reactive physical processes, which is otherwise considered computationally impossible using the conventional DSMC method. This conservative weighting scheme (CWS) improves greatly the statistical uncertainties by decreasing the weighting factors of trace-species particles and ensures the conservation of both momentum and energy between two colliding particles with large dierence of weighting factors. This CWS is further extended to treat reactive processes for gas-phase and surface reactions with trace species, which is called extended conservative weighting scheme (ECWS). A single-cell equilibrium simulation is performed for verifying both the CWS and ECWS in treating trace species. The results of using CWS show that it is most ecient and accurate for weight ratio (trace to non-trace) equal to or less than 0.01 for ows with two and three species. The results of a single-cell simulation using ECWS for gas-phase reaction and surface reactions show that only ECWS can produce acceptable results with reasonable computational time. Copyright ? 2003 John Wiley & Sons, Ltd.
KEY WORDS: LPCVD; conservative weighting scheme; extended conservative weighting scheme; trace species; reactive; single-cell simulation
1. INTRODUCTION 1.1. Fundamentals of CVD process
Chemical vapour deposition (CVD) is one of the principal unit operations in microelectronics processing [1]. It has been applied in many ways such as the formation of thin lm of metals, semiconductors and insulating materials of high purity, encompassing a wide range of material properties. The basic principle of CVD is that a suitable combination of reactant gases is
∗Correspondence to: J.-S. Wu, Department of Mechanical Engineering, National Chiao-Tung University, 1001
Ta-Hsueh Road, 30050 Hsinchu, Taiwan.
†E-mail: chongsin@cc.nctu.edu.tw ‡Associate professor.
brought in contact with a wafer surface that is maintained at an elevated temperature. Gas-phase precursor(s) may or may not react or decompose homogeneously into dierent gas-Gas-phase species depending upon the operating conditions. One or more of the gas-phase species reacts heterogeneously at the surface to deposit a solid lm. Usually the activation energy of the decomposition process is provided by the thermal energy of the wafer. It is generally agreed that a number of required steps, which must take place in the deposition sequence aecting the lm growth [2], including: (a) Transport of precursor(s) to the deposition zone; (b) Surface adsorption of the precursor(s); (c) Transformation of the precursor(s) to lm constituent; (d) Incorporation of lm constituent into crystal lattice (lm growth); (e) De-sorption of reaction by-products; (f) Transport of by-products away from the deposition zone.
1.2. Modelling of CVD process
The studies about modelling CVD process can be divided into two major categories. One is the study of the macro-scale analysis. This is about the whole ow led in the reactor, i.e. the spatial distribution of velocity, energy and density from the viewpoint of reactor size. These were evaluated by the Navier–Stokes equations [3–7] or the DSMC method [8–11]. The use of Navier–Stokes equations is most eective when the ow conditions are well in the hydro-dynamic regime. In contrast, the DSMC method [18] is ecient in the slip and transitional regimes, and is the only viable method in transitional regime. In addition, it is much easier to incorporate micro-scale physical models in the DSMC method. If the ow conditions are in the transitional regime, DSMC is the most popular method employed in this category, e.g. References [12, 13]. However, in the actual CVD process these physical phenomena strongly interact each other and the multi-scale analysis must be needed. For example, the angular distribution of the precursor molecular velocity and energy at the wafer surface is critical in determining the step coverage quality in trench or via holes. This information can be provided through the reactor-scale simulation. In addition, gas-phase chemistry inuences the precursor compositions at the wafer surface, which in turn changes the uniformity of deposited lm at the wafer surface. Therefore, a numerical approach suited to treat multi-scale process in CVD is required. Indeed, the direct simulation Monte Carlo (DSMC) method is one of the best choices among the candidates.
DSMC is a physically accurate particle method for the computation of non-equilibrium gas ows. The technique is most useful in circumstances where there are insucient numbers of collisions in the ow to maintain the equilibrium forms of the distribution functions describ-ing the various energy modes of the gas. Generally, such conditions prevail when the average distance between successive collisions of each particle, the mean free path, is comparable to the characteristic length scale of the ow. This type of non-equilibrium condition occurs in a variety of problems of current interest. These include hypersonic ows around vehicles ying at high altitude in planetary atmospheres, e.g. Reference [14], ows from small rockets used on satellites for control, e.g. Reference [15], ows involved in the synthesis of thin lms, e.g. Reference [16], and ows in micro-scale mechanical structures, e.g. Reference [17], to name a few.
Since the particle-based nature of the DSMC method, it is well suited to handle the ows and reactions in CVD process. Applying the DSMC method to simulate CVD process is generally divided into two parts, i.e. gas-phase (homogeneous) and surface (heterogeneous) reactions. In the gas-phase reaction, collisions are handled on a probabilistic basis. Within
each cell, potential collision partners are randomly selected without regard to their relative positions. A collision is performed if their collision probability exceeds some random fraction. And so does the chemical reaction if the steric factor exceeds another random fraction [18]. This process is repeated until the correct collision rate is obtained in each cell. Only binary collisions are considered as the diluteness of the gas makes any three-body collisions highly improbable. Additionally, in the surface reaction, the reaction can be simulated depending on the reactive probability,, using the same concept as in the gas-phase. This reactive probability may depend on the species concentration near the wafer surface [3, 5] or impinging energy of the precursor [19].
It is unfortunate that the chemical species of most importance in CVD process often occur in very small quantities. The amount of the trace gas, for example, SiH2, existing in CVD is usually very small for low-pressure environment, on the order of 10−4 or even smaller for the mole fraction. Similar situation occurs in plasma-enhanced CVD ows. This presents a major diculty to the DSMC algorithm. To simulate one particle of the trace gas in a computational cell at this mole fraction would require simulation of 10 000 other particles. This approach would therefore require at least 10 million particles or even more for a simu-lation of modest size. Even with a numerically ecient DSMC code implemented on parallel computers, this simulation would require hundreds of hours or even more of execution time. To circumvent this diculty, a conservative weighting scheme (CWS) was developed by Boyd [20] for the DSMC technique in which the physical weight of a particle W depends upon its chemical species. Thus, a particle representing a trace species would be given a lower weight than particles representing more abundant species. The scheme can also be extended to treat the chemical reaction for the steric factor in gas-phase and the reactive probability of the surface reaction, which we call extended conservative weighting scheme (ECWS) [21].
In summary, the objectives of the current research are listed as follows: (1) to implement and to assess a conservative weighting scheme for the DSMC method treating the trace species non-reactive gas ow; (2) to implement and to assess an extended conservative weighting scheme to deal with both gas-phase and surface reactions in a model CVD process.
2. NUMERICAL METHOD 2.1. The DSMC method
The basic idea of DSMC is to calculate practical gas ows through the use of no more than the collision mechanics. The molecules move in the simulated physical domain so that the physical time is a parameter in the simulation and all ows are computed as unsteady ows. An important feature of DSMC is that the molecular movement and intermolecular collisions are uncoupled over the time intervals that are much smaller than the mean collision time. Both the collision between molecules and the interaction between molecules and solid boundaries are computed on a probabilistic basis and, hence, this method makes extensive use of random numbers. In most practical applications, the number of simulated molecules is extremely small compared with the number of real molecules. In general, the procedure of DSMC method consists of four major steps: moving, indexing, collision and sampling. In the current study, the variable hard sphere (VHS) model [18] and no time counter (NTC) [18] are used to
simulate the molecular collision kinetics. The details of the procedures and the consequences of the computational approximations are well documented in Bird’s monograph [18] and thus are not repeated here.
2.2. Conservative weighting scheme (CWS) for non-reactive ows
Under equilibrium conditions, a species-dependent weighting scheme proposed by Bird [22] will approximately conserve linear momentum and energy over a suciently large number of collisions. However, the scheme does not conserve these properties explicitly at each collision as clearly demonstrated by Boyd [20]. Indeed, Bird [22] recommended against the use of this scheme. Hence, there is a requirement for developing a weighting scheme
for the DSMC technique that does conserve explicitly linear momentum and
energy.
To overcome the problem mentioned above, Boyd [20] proposed a conservative weighting scheme, which is described briey in the following. The rst stage of the conservative weight-ing scheme is to split the particle of abundant species (W1) into a particle with weight W2 (trace species) and a particle with weight of W1–W2 when two particles (trace and abundant species) collide. Then, a collision is performed using the conventional DSMC procedure for the two particles that have the same weight W2. The nal stage is to merge together the two particles that were split such that the each linear momentum in three physical directions is exactly conserved. Unfortunately, it does not explicitly conserve total energy. But the energy dierence (loss) caused by this split-merge process is found to be proportional to the weight ratio W2=W1 (¡1). Thus, the conservative weighting scheme proposed by Boyd [20] nearly conserves total energy as this weight ratio is much smaller than unity. The split-merge process described in the above can be summarized as Figure 1.
It was argued that if the split-merge scheme is employed at each collision, then energy is continuously lost from the system because of energy loss [20]. Boyd [20] also proposed some practical remedies to keep this energy loss to a minimum by adding lost energy to the central-mass energy in a subsequent collision. In general, this energy should only be added to collisions between particles having the maximum weight used in the simulation to keep this eect a minimum; that is, between two non-trace particles (leading to a non-trace collision). Thus, energy conservation is essentially maintained for each iterative step of the simulation.
2.3. ECWS for chemical reactive ows with trace species
In the following, we will briey describe the homogeneous gas-phase and heterogeneous surface chemical reaction involved in the CVD process (e.g. the silicon vapour deposition) and the details of corresponding ECWS we will use to handle these processes.
2.3.1. ECWS for homogeneous gas-phase reaction. Consider a model gas-phase reaction as follows:
AB + M↔A + B + M (1)
where AB is the trace species (e.g. SiH4), M is any abundant species (e.g. H2), A is the most trace species (e.g. SiH2) and B is the decomposed gas (e.g. H2) in the current study.
W2m1 W2m2 W1m1 (W1-W2)m1 non-collision collision W2m1 W2m2 W2m1 W1m1 (W1-W2)m1 Splitting Step Colliding Step Merging Step (W1-W2)m1
Figure 1. Schematic diagram of CWS for non-reactive ow.
In the following, we will describe how to implement ECWS considering chemical reactions in DSMC simulation. The basic idea for treating homogeneous decomposition is to decrease gradually the weights of reactants for every reactive collision [2]. Figure 2 shows the schematic diagram of ECWS for treating trace species in reactive ows. In summary, ECWS consists of three main steps as described in the following.
Splitting step: The rst step is to divide AB particle and M particle into two parts, re-spectively, as shown in Figure 2. The one with superscript ‘r’ is the part of reactive colli-sion, while the other with superscript ‘n’ represents the part of non-reactive collision. If we
Wr
MmM WrABmAB
WnMmM W n
ABmAB
reactive collision non-reactive collision
Wr ABmAB WnMmM WnABmAB WAmA Wr MmM Wr MmM WBmB WnMmM WnABmAB Wr MmM W n MmM WMmM M AB M Splitting Step Colliding Step Merging Step WMmM WABmAB
Figure 2. Schematic diagram of Extended CWS for treating homogeneous decomposition with trace species.
assume that a particle AB is a collection of NAB (=WAB=WA) equivalent particles (one equiv-alent particle is equal to one particle A), then the number of equivalent particles involved in each chemical reaction is
Nr
AB=WAB=WAPr=NABPr (2)
where Pr is the reaction probability, where the format depends upon the model used (e.g. TCE model [18] as used in the current study unless otherwise specied). This means that Nr
AB equivalent particles out of NAB particles are involved in each reaction. Furthermore, the weighting factor for reactive and non-reactive part of AB (or M) particle can be expressed
respectively as
Wr
AB=WMr=NABr (3)
Wn
AB=WAB−NABr; WMn=WM−NABr (4)
Note that mass conservation is strictly satised during the splitting process. In addition, Wr
AB (=WMr), WABn and WMn vary depending upon the ow conditions during each collision. Colliding step: Collisions are then performed using the conventional DSMC method for both parts as illustrated in the middle part of Figure 2. First, we consider the reactive part of the collision. The reaction model [e.g. TCE model in the current study] can be applied with ABr and Mr having the same weight, and species A and B are created with the weight factor WA andWB, respectively. The number of the created A and B particles is WABr . Though equivalent particles ABr vanishes, the other equivalent particleMr still remains with the weight factor of Wr
M. Non-reactive part of the collision between particles with weight factors of WMn and WABn is handled using CWS, which is the same as that described in previous section. Therefore, all the post-collision states (velocities and energies) of all particles can be determined. In the current method, mass conservation is, however, enforced by decreasing the weight factor of AB particle, rather than by consuming the number of AB particles.
Merging step: The last step is to determine the velocity and energy of third body M as shown in Figure 2. Two equivalent particles, Mr andMn, have the velocity vectors, CMr and Cn
M, having the dierent weights,WMr and WMn, respectively. Again, by applying CWS, we can merge the two equivalent particles. The nal velocity vector CM is evaluated as
C
M= (1−)CMn +CMr (5)
where = NABr =WM. The energy loss between the pre- and post-collision states can be readily derived as
E =
i mM(1−)C n
iCir (6)
This energy loss is added by the same way as CWS. It is noted that this energy loss is small if is much less than unity, which is always the case in the reaction of silane decomposition. The inelastic collision involving internal energy can be determined similarly.
It is noted that each particle of products A and B will be created each time step with the appropriate weight factors and the conservation of momentum (exactly) and energy (nearly) is satised at each collision by ECWS.
The conventional formulation of the steric factor (TCE model) for the VHS gas model rst proposed by Bird [23] depends on the internal contributed degree of freedom as
Pr Ec kT = 1=2T ref 2ref(kTref)−1=2+w ( + 2−w) ( + + 3=2) mr 2kTref 1=2(E c−Ea)++1=2 E+1−w c (7) where has a value of one for like and two for unlike collision partners, and are constants,ref denes a reference cross section at a temperatureTref,k is Boltzmann’s constant, w = !ab−1=2, where !ab is the viscosity–temperature exponent that is specied for collisions between molecules of species A with those of species B, is the average number of internal degrees of freedom, and Ea is the activation energy of the reaction. For procedures of gas-phase particle collisions considering chemical reactions, see Figure 3 for details.
No Yes
No Yes
Calculate possible collision pair numbers
Select collision pairs and determine if collision occurs
Handle the chemistry reaction by split-merge method
Assign energy and velocity to the product molecule and remove the
reactive molecule Determine if redistribution of V-T
or R-T energy occurs respectively
Assign energy according to L-B
distribution
Assign molecules velocities
Finish all possible collisions
Next time step If ETc > Ea
Yes No
If chemistry occur
Figure 3. Sketch of ow chart for gas-phase particle collision.
2.3.2. ECWS for heterogeneous surface reaction. Consider two model surface reactions:
AB→C + 2B (8)
A→C + B (9)
where, for example, A stands for SiH2, B for H2, AB for SiH4 and C for Si in silane-based silicon vapour deposition. To simplify the calculation, we consider the model of Moat [6], a non-linear Arrhenius temperature dependence, for the surface reaction probability of AB (e.g. silane). In addition, we set the surface reaction probability for the highly reactive open shelled molecule A (e.g. SiH2) to one as presented in Reference [6].
In the DSMC simulation, if is much smaller than unity, enormous computational time is required to achieve reasonable sample size on the surface. Similar to the gas-phase reaction, we use extended CWS to circumvent these problems. Consider the case that a computational
A Step 1 WAmA A* Step 2 A Escaped (1- )WAmA A WAmA Adsorbed WAmA
Figure 4. Schematic diagram of Extended CWS for treating heterogeneous reaction.
species A with weight WA being adsorbed to the surface with the probability . If we assume that a computational particle A is a set of WA equivalent particles, the number of equivalent particles that are adsorbed to surface every time the particle A hits the surface is expected as
WS
A =WA (10)
where the superscript S represents the surface related properties. Then the computational particle A is divided into two parts again, as shown in Figure 4. That is, one part is adsorbed to the surface with weight WS
A , and another escapes from the surface with weight WA–WAS which enforces the conservation of mass.
3. RESULTS AND DISCUSSIONS
In the following, we have performed the single-cell DSMC simulation for non-reactive and reactive ows by using CWS, variable weighting scheme by Bird [22] and constant weighting scheme, respectively. By comparing the velocity distribution, chemical reaction rate coe-cient and deposition distribution, it can be clearly shown that CWS is much superior to the conventional scheme if trace species is involved.
3.1. Simulation for two-component mixture without chemical reactions
The rst simulation includes two species: Ar and He with dierent weight ratios. In this sim-ulation, we would like to know the eectiveness and accuracy of the conservative weighting scheme for non-reactive ows. We perform the simulation with dierent weight ratios using variable weighting scheme by Bird [22], constant weighting scheme and CWS, respectively. The total simulation particles are 10 000 in the cell, and the weight ratios are WAr=WHe= 0:1 (10%), 0.05 (5%) and 0.01 (1%), respectively. Note that the number in the parenthesis repre-sents mole fraction of Ar. Thus, Ar is a potential trace species in this case. The temperature in the cell is preset as 1000K. Figures 5 and 6, with weight ratios of 1=9 and 1=19, respectively, show the velocity distribution after 20 000 time steps. It is clear that the conventional scheme (constant weighting) performs better than both the variable weighting scheme by Bird [22] and CWS. This is because CWS has energy loss from non-trace gas, and we must add the lost energy to non-trace gas. When the weight ratio is not small enough, the energy loss will be relatively large, and hence it will aect the energy distribution. In turn, the velocity distribution will deviate strongly from the expected M–B (Maxwell–Boltzmann) distribution, which it should be in the current test case. But when the weight ratio is equal to or less than 0.01, as shown in Figure 7, the CWS can perform as accurately as the conventional scheme (constant weighting) once the weight ratio is equal to or less than 0.01, while the variable weighting scheme by Bird [22] shows some over-sampling at the velocity of zero. Figures 8–10 illustrate the relative error as a function of simulation time for the two-species case withWAr=WHe= 0:1, 0.05 and 0.01, respectively. Note that the relative error is dened as the root mean square of the sum of deviation from the M–B distribution relative to the M–B
-4000 -2000 0 2000 4000 Velocity (m/s) 0.00 0.20 0.40 0.60 0.80 1.00 Distribution
SYM. SPECIES* SCHEME
Ar Conventional He Conventional Ar CWS He CWS Ar WS (Bird) He WS (Bird) Ar M-B Distr. He M-B Distr. *n(Ar) : n(He) = 1 : 9
-4000 -2000 0 2000 4000 Velocity (m/s) 0.00 0.10 0.20 0.30 0.40 Distribution
SYM. SPECIES* SCHEME
Ar Conventional He Conventional Ar CWS He CWS Ar WS (Bird) He WS (Bird) Ar M-B Distr. He M-B Distr. *n(Ar) : n(He) = 1 : 19
Figure 6. Velocity distribution with weight ratio 1:19.
-4000 -2000 0 2000 4000 Velocity 0.00 0.04 0.08 0.12 0.16 0.20 Distribution
SYM. SPECIES* SCHEME
Ar Conventional He Conventional Ar CWS He CWS Ar WS (Bird) He WS (Bird) Ar M-B Distr. He M-B Distr. *n(Ar) : n(He) = 1 : 99
0 4000 8000 12000 16000 20000 Time Step 0.00 1.00 2.00 3.00 4.00 5.00 6.00 Relative Error
SYM. SPECIES* SCHEME Ar Conventional He Conventional Ar CWS He CWS Ar WS (Bird) He WS (Bird) *n(Ar) : n(He) = 1 : 9
Figure 8. Relative error as a function of the number of simulation time steps with weight ratio 1:9.
0 4000 8000 12000 16000 20000 Time Step -0.40 0.00 0.40 0.80 1.20 1.60 Relative Error
SYM. SPECIES* SCHEME Ar Conventional He Conventional Ar CWS He CWS Ar WS (Bird) He WS (Bird) *n(Ar) : n(He) = 1 : 19
0 4000 8000 12000 16000 20000 Time Step 0.00 0.10 0.20 0.30 0.40 Relative Error
SYM. SPECIES* SCHEME
Ar Conventional He Conventional Ar CWS He CWS Ar WS (Bird) He WS (Bird) *n(Ar) : n(He) = 1 : 99
Figure 10. Relative error as a function of the number of the simulation time steps with weight ratio 1:99.
distributions by dividing the velocity range (−5000 –5000m=s) into 50 sections. In the case of weight ratio 0.01, the relative errors of both species by the CWS decrease to a reasonably low value in a few time steps, while that of trace species by the conventional method maintains unacceptably large as simulation continues. In the long run, the variable weighting scheme [22] achieves about the same level of errors as the CWS, however, it progresses with much higher errors in the early stage of simulation, which should be caused by the non-conservation of momentum and energy in each collision mentioned previously. This result shows not only the CWS is as accurate as the constant weighting scheme but is also very eective for dealing the ow involving trace species with the weight ratio equal to or lower than 0.01.
3.2. Simulation for three-component mixture without chemical reactions
This simulation includes three species, Ne, Ar and He, with weight ratio 1:99:9900 in order. The simulation conditions are similar to those of the two-species case. The total simulation particles increase up 100 000, and the temperature in the cell is kept at 1000K. Figure 11 shows the velocity distribution after 20 000 time steps for all three schemes. It clearly demonstrates that the conventional scheme (constant weighting) is incapable of handling the trace species in the period of simulation time because very few collisions of Ne–Ne occur in the cell. Again the velocity distribution by the variable weighting scheme [22] approximately ts the M–B curve, but with an over-sampling at the velocity of zero, which is similar to that in two-species simulation. Figure 12 shows the relative errors of both the trace species (Ar) and the most trace species (Ne) by the conventional scheme (constant weighting) are too large and decrease very slowly. As for the variable weighting scheme [22], similar situation as that
-4000 -2000 0 2000 4000 Velocity (m/s) 0.00 0.04 0.08 0.12 0.16 0.20 Distribution
SYM. SPECIES* SCHEME
Ne Conventional Ar Conventional Conventional He Ne CWS CWS Ar He CWS Ne WS (Bird) WS (Bird) WS (Bird) Ar He Ne M-B Distr. M-B Distr. M-B Distr. Ar . He *n(Ne):n(Ar):n(He) = 1:99:9900
Figure 11. Velocity distribution with weight ratio 1:99:9900.
0 400 800 1200 1600 2000 Time Step 0.00 0.20 0.40 0.60 Relative Error
SYM. SPECIES* SCHEME
Ne Conventional Ar Conventional He Conventional Ne CWS Ar CWS He CWS Ne WS (Bird) Ar WS (Bird) He WS (Bird) *n(Ne):n(Ar):n(He) = 1:99:9900
Figure 12. Relative error as a function of the number of the simulation time steps with weight ratio 1:99:9900.
of two-species simulation occurs. Thus, it is clear that the CWS performs much better than the conventional schemes if trace species are involved in DSMC simulation.
3.3. Simulation with gas-phase reaction
In this simulation, we would like to know if ECWS could reproduce the correct chemical reaction rate for before its application to the modelling of practical reactive ows. The gas is a mixture of SiH4 (1%, trace species) and H2 (99%) initially in a single cell. The initial total simulation particles are 20 000. The temperature in the cell varies from 800 to 1500K, which will change the reaction rate accordingly. We use total collision energy model (TCE) proposed by Bird [18]. Note that the variable weighting scheme [22] is not performed in the following simulations for reactive ows. It lies in the fact that the large number of simulated molecules at the same position and with the same velocities (once reaction (1) occurs in the variable weighting scheme) will denitely cause serious bias problems in the collisional phase in a statistical simulation like DSMC, in addition to the non-conservation of momentum and energy during collision.
Figure 13 shows the chemical reaction rate coecient Kf calculated by TCE model using theoretical Arrhenius law (related coecients from Kleijn [3]), conventional scheme (constant weighting) and CWS as a function of temperature. Note that total internal degree of freedom (IDOF) is selected as 11. Both the conventional scheme and ECWS agree well with the theoretical value in the range of 1000 –1500K. It is, however, clear that the number of time steps to accumulate reasonable sample size for the conventional scheme is much larger than that of CWS (about two order of magnitudes larger), as shown in Figure 14. It becomes nearly impossible to use the conventional scheme for practical computations especially at lower temperatures, (800 and 900K). In the current case, data of the conventional scheme are not shown (Figures 13 and 14) due to the inhibitive long computational time.
Figure 15 shows the variations of the predicted rate coecient as a function of simulation time steps at 1200K (IDOF= 11). It is clear that the rate coecient computed by the con-ventional scheme does not coincide with theoretical value even after a long period of time,
800 900 1000 1100 1200 1300 1400 1500 Temperature (K) 1e-026 1e-025 1e-024 1e-023 1e-022 1e-021 1e-020 Rate Coefficient Kf SYM. SCHEME Theory Conventional ECWS
800 900 1000 1100 1200 1300 1400 1500 Temperature (K) 1000 1e+004 1e+005 1e+006 1e+007 1e+008 Time Step SYM. SCHEME Conventional ECWS
Figure 14. Number of time steps required to approximate theoretical value of rate coecient (IDOF= 11).
0 1000000 2000000 3000000 4000000 Time Step 0 1e-022 2e-022 3e-022 4e-022 5e-022 Rate Coefficient Kf SYM. SCHEME Theory Conventional Scheme ECWS
Figure 15. Rate coecient Kf as a function of simulation time steps calculated (conventional scheme and ECWS at 1200K,IDOF= 11).
Table I. Rate coecient Kf computed by CS and ECWS, and simulation time steps and associated
errors as compared with theoretical value.
Temp. Kf Kf Error (%) Time steps Kf Error (%) Time steps
(K) (Theory) (CS∗) (CS) (CS) (ECWS†) (ECWS) (ECWS)
800 2:40E − 26 — — — 2:17E − 26 −9:58 6 000 000
900 4:91E − 25 — — — 4:15E − 25 −15:48 1 900 000
1000 5:49E − 24 3:03E − 24 −44:81 20 000 000 5:68E − 24 3:46 160 000
1100 3:96E − 23 4:01E − 23 1:26 5 600 000 4:29E − 23 8:33 60 000
1200 2:05E − 22 1:62E − 22 −20:98 1 400 000 2:04E − 22 −0:49 55 000
1300 8:26E − 22 9:69E − 22 17:31 600 000 8:21E − 22 −0:61 27 000
1400 2:73E − 21 2:69E − 21 −1:47 300 000 2:64E − 21 −3:30 10 000
1500 7:61E − 21 7:90E − 21 3:81 210 000 7:60E − 21 −0:13 6000
∗CS = conventional scheme.
†ECWS = extended conservative weighting scheme.
1000 1100 1200 1300 1400 1500 Temperature (K) 1e-024 1e-023 1e-022 1e-021 1e-020 Rate Coefficient Kf
ECWS (Internal DOF = 11)
Theory TCE
Boyd & Stark (1990)
Figure 16. Rate coecient calculated by ECWS, internalDOF= 11.
e.g. 4 million time steps due to the very low collision probability in the current test case. The rate coecient, predicted by ECWS, however, reaches the theoretical value very fast and stays at the approximately constant value afterwards as simulation continues. Table I presents the errors and time steps required for sampling by the conventional scheme and ECWS. It can be seen clearly that ECWS is superior to the conventional scheme for simulating chemical reaction involving trace species.
To identify the eects of IDOF using ECWS, Figures 16 and 17 present the rate coecients (by ECWS) as a function of the temperature for total internal degree of freedom of 11 and 3,
1000 1100 1200 1300 1400 1500 Temperature (K) 1e-024 1e-023 1e-022 1e-021 1e-020 Rate Coefficient Kf
ECWS (Internal DOF = 3)
Theory TCE
Boyd & Stark (1990)
Figure 17. Rate coecient calculated by ECWS, internalDOF= 3.
respectively. The results show that the prediction of rate coecient by TCE along with the model proposed by Boyd and Stark [19] for IDOF= 11 are both in good agreement with the theoretical value in the temperature range of 1000 –1500K. However, it deviates from the theoretical value if IDOF= 3 is used, especially in the higher temperature range (1400 – 1500K). This may be caused by the increasing importance of internal degree of freedom at higher temperature. This has also been observed and conrmed by Hsu [24], where he has shown that DSMC simulation, using internalDOF= 11, can reproduce the theoretical rate coecient in the similar temperature range.
3.4. Simulation with surface reaction
The number of initial simulation particles in a single cell is 100 000, and the temperature is 1500K in the gas as well as at wall. Only the bottom wall is able to have surface reactions. It is divided into 25 sections for counting the deposited particles. Figure 18 shows the total number of deposited particles on the surface at dierent surface reaction probability after 20 time steps of simulation. It is clear that ECWS produces about 100 times as many deposited particles compared with the conventional scheme within 20 time steps of simulation for all the considered. This improvement reduces greatly the statistical uncertainties of evaluating deposition uniformities on the surface otherwise caused by the low surface reaction probability.
Figure 19 presents the number of simulation time steps required as a function of the surface reaction probability within a preset uniformity tolerance. This tolerance is dened as the sum
1E-006 1E-005 1E-004 1E-003 1E-002
Surface Reaction Probability γ
0.1 1 10 100 1000 10000 100000 Deposition Particles SYM. SCHEME Conventional ECWS
Figure 18. Number of deposited particles on the surface as a function of surface reaction probability , after 20 time steps.
1E-006 1E-005 1E-004 1E-003 1E-002
Surface Reaction Probability γ
1 10 100 1000 10000 100000 1000000 10000000 Time Steps SYM. SCHEME Conventional ECWS
Figure 19. Number of time steps to obtain acceptable surface reaction probability (normalized error ¡0:8).
of the normalized absolute dierences (between average and local value) for all deposited sites on the surface. In the simulation, the preset tolerance is set as 0.8. It is clear that the conventional scheme generally requires two order-of-magnitude large number of time step to achieve preset uniformity on the surface.
-0.4 -0.2 0 0.2 0.4 Position (m) 0 0.4 0.8 1.2 1.6 2 2.4
Normalized Deposited Particle Number
SYM. SCHEME
Conventional ECWS
Figure 20. Spatial distribution of normalized deposited particle number for = 10−3.
Figure 20 illustrates the spatial distribution of normalized deposited particle number along the surface. It is again clear that the uniform distribution of deposited particle number can be more or less reproduced using ECWS. On the contrary, the results become unacceptable using the conventional scheme for case of low surface reaction probability.
In summary, the above verications have clearly demonstrated the superiority of CWS over the conventional scheme in simulating chemical reaction process (both homogeneous and heterogeneous) involving trace species.
4. CONCLUSIONS
In the current study, the conservative weighting schemes (CWS and ECWS) are assessed to handle the trace species in the DSMC simulation with and without chemical reaction using a single-cell simulation. In summary, the major ndings of the current research are listed as follows:
1. CWS is shown to be eective and accurate in DSMC simulation without chemical reac-tion for two-species and three-species mixtures with concentrareac-tion ratio (trace to abun-dant) equal to or less than 0.01, which the conventional scheme (constant weighting) is dicult to use practically and the variable weighting scheme by Bird [22] is comparably inaccurate in short period of running time.
2. ECWS is shown to be eective and accurate in DSMC simulation with gas-phase reac-tions with trace species and surface reacreac-tions with very small reaction probability, which are both otherwise impossible to model using conventional scheme.
3. ECWS using TCE model and that proposed by Boyd and Stark [19] in DSMC simulation produce similar results in simulating the chemical reaction rate for silane decomposition in gas phase.
4. In the temperature range of 800 –1500K, satised results of rate coecients by DSMC simulation are obtained if IDOF= 11 is used for silane decomposition.
From the results of the current study, it can be concluded that the ECWS is highly potential in simulating reactive ows using DSMC involving gas-phase if there is any critical trace species present and surface reactions if low surface reaction probability involved. Note that the ECWS can be applied to treat plasma-enhanced CVD ows involving critical trace species in a similar way. In addition, simulating the chemical deposition onto a single wafer by a downward impinging jet (e.g. a typical LPCVD chamber), using DSMC incorporating ECWS, is currently in progress and will be reported in the near future.
REFERENCES
1. Hitchman E, Jensen K. Chemical Vapour Deposition: An Overview in Chemical Vapour Deposition: Principles and Applications. Hitchman M, Jensen K (eds). Academic Press: New York, 1993.
2. Coronell DG. Simulation and analysis of rareed gas ows in chemical vapour deposition processes. Ph.D. Thesis, MIT, 1993.
3. Kleijn CR. A mathematical model of the hydrodynamics and gas-phase reactions in silicon LPCVD in a single-wafer reactor. Journal of the Electrochemical Society 1991; 138(7):2190–2200.
4. Coltrin ME, Kee RJ. A mathematical model of the gas-phase and surface chemistry in GaAs MOCVD. Material Research Society Symposium Proceedings 1989;145:119.
5. Coltrin ME, Kee RJ, Evans GH. A mathematical model of the uid mechanics and gas-phase chemistry in a rotating-disk chemical vapour deposition reactor. Journal of the Electrochemical Society 1989;136(3):819–829. 6. Moat HK, Jensen KF. Three-dimensional ow eects in silicon CVD in horizontal reactors. Journal of the
Electrochemical Society 1988;135(2):459–471.
7. Wang YB, Chaussavoine C, Teyssandier F. 2D modelling of a non-conned circular impinging jet reactor; Si chemical vapour deposition. Journal of Crystal Growth 1993;126:373–395.
8. Coronell DG, Jensen KF. Simulation of rareed gas transport and prole evolution in nonplanar substrate chemical vapour deposition. Journal of the Electrochemical Society 1994; 141(9):2545–2551.
9. Nanbu K, Igarashi S, Watanabe Y. Molecular simulation of lm growth rate in the low-pressure CVD method. JSME International Journal Series B-Fluids & Thermal Engineering 1990;56(524):892–897.
10. Nanbu K, Mitamura S, Igarashi S. Growth rate of silicon lms on a wedge-shaped substrate in the LP-CVD reactor. JSME International Journal Series B-Fluids & Thermal Engineering 1991;57(542):3526–3530. 11. Nanbu K, Sugawara T, Igarashi S. Monte Carlo simulation of the growth rate of lms in a CVD diusion
reactor. JSME International Journal Series B-Fluids & Thermal Engineering 1991;57(543):3760–3764. 12. Ikegawa M, Kobayashi J. Semiconductor deposition prole simulation using direct simulation Monte Carlo
method. JSME International Journal Series B-Fluids & Thermal Engineering 1993;59(567):3365–3371. 13. Masato I, Junichi K. Deposition prole simulation using the direct simulation. Journal of the Electrochemical
Society 1989; 136(10):2982–2986.
14. Nance RP. Monte Carlo simulation of three-dimensional hypersonic ows on parallel architectures. Master Thesis, North Carolina State University, 1995.
15. Kannenberg KC. Computational method for the direct simulation Monte Carlo technique with application to plume impingement. Ph.D. Thesis, Cornell University, Ithaca, NY 1998.
16. Plimpton S, Bartel T. Parallel particle simulation of low-density uid ows. U.S. Department of Energy Report No. DE94-007858, 1993.
17. Wu JS, Tseng KC. Analysis of micro-scale gas ows with pressure boundaries using direct simulation Monte Carlo method. Computers & Fluids 2001;30:711–735.
18. Bird GA. Molecular Gas Dynamics and the Direct Simulation of Gas Flows. Oxford University Press: New York, UK, 1994.
20. Boyd ID. Conservative species weighting scheme for the direct simulation Monte Carlo method. Journal of Thermophysics and Heat Transfer 1996;10(4):579–585.
21. Sakiyama Y, Takagi S, Matsumoto Y. Full simulation of silicon chemical vapour deposition process. In 22nd International Symposium on Rareed Gas Dynamics, Bartel TJ, Gallis MA (ed.). 2000; 206 –213.
22. Bird GA. Molecular Gas Dynamics. Oxford University Press: Oxford, England, UK, 1976.
23. Bird GA. Monte Carlo simulation in an engineering context. In Progress in Astronautics and Aeronautics: Rareed Gas Dynamics, vol. 74, Pt. 1, Sam S. Fisher (ed.). AIAA: New York, 1981; 239 – 255.
24. Hsu CC. Molecular simulation of monosilane gas-phase decomposition and diusion by DSMC. Master Thesis, Department of Chemical Engineering, National Taiwan University of Science and Technology, 1997.
![Figure 13 shows the chemical reaction rate coecient K f calculated by TCE model using theoretical Arrhenius law (related coecients from Kleijn [3]), conventional scheme (constant weighting) and CWS as a function of temperature](https://thumb-ap.123doks.com/thumbv2/9libinfo/7507211.117181/15.850.195.650.684.930/reaction-calculated-theoretical-arrhenius-conventional-weighting-function-temperature.webp)
