A plasma chemistry model for H2/SiH4 mixtures used in PECVD processes

Plasma chemical processes in H2/SiH4 discharges are critically reviewed. A model set of reactions is proposed which includes temperature and pressure-dependent reaction rates and describes SiyHx (y ≤ 3) chemistry. Using a 2D fluid plasma simulator, the model has been tested under three different set of operating conditions. First, it has been validated against the experimental benchmark data (Horvath and Gallagher (2009) J. Appl. Phys. 105, 13304). Based on considerations of atomic hydrogen content, the branching of SiH4 dissociation channels and the H surface loss probability have been defined more accurately. Then, simulations have been also performed for the plasma source of a PECVD tool from Meyer Burger Germany. A very good agreement between the computed and experimentally determined deposition rates can be stated.


Introduction
Thin films of hydrogenated amorphous silicon a-Si:H are widely used in photovoltaic industry to produce efficient solar cells [1]. In order to become more competitive, the industrial PECVD processes producing these layers require optimization in terms of film properties, uniformity, high throughput.
Detailed understanding of the physics and chemistry behind the deposition process is indispensable in realizing these goals. Comprehensive modelling and simulation of a deposition device can be extremely helpful in achieving such understanding.
A lot has been already done in the research field of SiH 4 /H 2 discharges. Back in the end of the 80s, the simulation studies by M. Kushner [2,3] laid the foundation in terms of full-scale modelling. The works by Perrin and co-authors [4][5][6][7][8][9] advanced the field further, unraveling various aspects of gas-phase and surface processes and providing data needed for modelling purposes. The summarizing work [7] offers compilation and critical review of the knowledge available in the mid-90s. After that, research in the field of a-Si:H deposition continued, but the focus shifted to plasma-surface interaction due to the rise of extensive computing. Such methods like molecular dynamics and kinetic Monte Carlo proved very useful in studying the details of the film growing mechanisms and radical-surface interactions, just to cite a couple of examples [10,11].
Despite all the progress, still lacking in the literature is a conclusive plasma chemistry model with consistent pressure-and temperature-dependent rate constants and an appropriate description of Si 2 H x radicals.
The present study seeks to fill the gap. The focus is on thorough compilation of reaction set and on updating the rate constants by making use of the literature which is beyond the scope of [7].
The paper is organized as follows. Following Introduction section, Plasma Chemistry Model presents full lists of species, reactions, and constants. Some crucial reaction pathways are critically reviewed. The described set of reactions has been used within the the fluid simulator Quantemol-VT, that is also briefly introduced. In Experiment section, we describe a-Si:H deposition performed by a capacitively coupled radio frequency plasma source in an industrial-type inline PECVD tool (Meyer Burger MAiA). This section also gives the design and operation of the TIMS experiment [12] whose data are used to validate our model. In Results and Discussions, a validation of our model is performed for three different cases of discharge conditions. The simulated and

Introduction
The goal of this work is to develop a model for H 2 /SiH 4 discharges which is appropriate in the range of operating conditions 0.1 Torr < p < 1.5 Torr and 300 K < T < 700 K. This means that the p and T dependencies of the reaction rates must be considered whenever possible. The film deposition under these conditions is essentially determined by the fluxes of neutral radicals of small Si-atom numbers [3]. Accordingly, our plasma model involves rather detailed chemistry of SiH x and Si 2 H x species. Si 3 H x molecules are represented by a single lumped state, which concept is briefly discussed in the end of this section. The relative contribution of ion fluxes is usually much smaller but may become more significant approaching the lower end of the pressure range [5]. Ion-neutral reactions would not generally change the total number of ions, so that there is no need in overly detailed volumetric ion chemistry. Therefore, only most principal ions are retained in the model, which are the primary products of the electron impact ionization of the precursor gases.
Species and reactions used in the model All the species and reactions between them are summarized in tables 1-3.
In total, there are 20 species in the model which are listed in table 1. The radicals and ions are supplied with respective surface loss probability ß and sticking coefficient s, following the notation from [7]. While ß describes the total losses on a boundary, sticking coefficient gives the portion of the incident flux which gets incorporated into the film, thus contributing to the film growth. The recombination coefficient γ represents the probability to form a stable volatile product, so that holds ß = s + γ. The choice of these probabilities is discussed in the end of this section.
The starting point of the plasma chemistry is the dissociation by electron impact, see reactions R1-R4 in table 2. For all the importance of the dissociation, there is no clarity on the branching of specific channels. There is common agreement that the main channels are R1 and R2, yet the ratio between them R1/R2 takes values in the wide range between 8.5 and 0.2 (e.g. [2]). This last extreme point R1/R2 = 0.17/0.83 is popular in the literature on modelling. It is based on a single experimental evidence, originating from the direct photolysis of silane by 147 nm radiation [30], and it was suggested afterwards that this might also relate to the dissociation by  [15], where the near-threshold behavior of the cross-sections have been carefully taken into account. This choice is also supported by our results discussed later. Note that, in   [15], the product hydrogen in reaction R2 is given as H 2 . This work remains an exception in this regard. We stick to 2H due to its wider acceptance in the literature (for example in [7,31]).

Pressure and temperature dependence of selected reactions
In the following, our choice of the rate constants of R20-R25 is explained in detail.
Magnitude at room temperature There have been numerous experimental determinations of the rate constant k 20 , all falling in the range between 2.1 × 10 −13 and 8.5 × 10 −12 cm 3 /s at room temperature. The most recent determinations (2.8 ± 0.3) × 10 −13 cm 3 /s [32] and (3.38 ± 0.16) × 10 −13 cm 3 /s [20] at 298 K agree well with each other. Moreover, in [20] was stated a good agreement with the value averaged across the literature (to the exclusion of evident outliers) (3.4 ± 0.5) × 10 −13 cm 3 /s, which seems to represent a well-established value at room temperature.

Temperature dependence
Being independent of pressure, the rate constant is sensitive to the temperature. This dependence was also measured in the two cited works. These data are in substantial agreement, showing only some divergence at the upper end of the temperature range. In view of this, the authors of [20] combine both datasets and provide a fit expression (1.37 ± 0.03) × 10 −10 exp[ -(1820 ± 10)/T] cm 3 /s for the temperature range 290-636 K. This temperature dependence is stronger than the expression cited in [7] 2.8 × 10 −10 exp[ −1250/T], giving discrepancies up to 40% in the range from 300 to 500 K. The latter expression is based on the results from [33], where the rate constant 4.3 × 10 −13 was only measured at room temperature, after which the Arrhenius parameters were proposed relying on some 'admittedly crude approximation' according to the authors' words. Besides, in [14] was argued that this rate constant in a real SiH 4 /H 2 plasma might be higher due to hot H atoms which are not fully thermalized after being produced by electron impact dissociation.

Magnitude at room temperature
The values of the rate constant found in the literature range from 2 × 10 −11 to 5 × 10 −10 cm 3 /s. In the analysis of a pulse radiolytic experiment [34], this constant was assumed 5 × 10 −10 as a first approximation. Further, it was found that, divided by two, this value gives better agreement with the experimental data. The data was acquired at 300 K and in the atmospheric pressure range. In [3] a value of 10 −10 was used. By modeling the results of their photolytic experiments [35], the authors inferred the rate constant (2 ± 1) × 10 −11 at 295 K and 9.5 Torr. Later, this value was quoted in [7]. However, in a later publication [8] with partly the same authors as in [7], a value of 5 × 10 −10 was involved in their chemistry model. Simulation studies of another photolysis experiment in [20] yielded a value of (4.3 ± 1.0) × 10 −11 cm 3 /s at 298 K and 500 Torr. Theoretical calculations [36] yielded 1.4 × 10 −10 cm 3 /s as association rate of SiH 3 + H at 300 K.
T and p dependence Ab initio theoretical calculations [37] demonstrated negligible pressure dependence (between 10 −3 and 100 bar at least) and only slight temperature dependence in the range 300-2000 K. The Arrhenius expression 6.9 × 10 11 T 0.736 exp(134.8/T) cm 3 /mol/s yields 1.2 × 10 −10 cm 3 /s at room temperature. Somewhat higher value of 2.1 × 10 −10 (T = 300 K) was obtained in another ab initio calculation study [21], that also confirmed missing pressure dependence in the range of interest and a weak temperature dependence. The interpolation expression for temperatures between 300 and 1000 K is 3.55 × 10 −10 T −0.05 exp(−77/T) cm 3 /s. We opted for the rate constant from the latter work because it also provides a consistent set of rates for H-abstraction (R22 and R23) from other SiH x radicals.
This pathway is of utter importance for SiH 2 radical and for production of higher silanes. It has two channels [23] : Si H R24 While channel (R24) is strongly dominating at lower temperatures, channel (R25) gains in importance at higher temperatures due to the opposite trend in T dependence. The rate constant k 25 (p,T) is given in figure 10(b) in [23]. At p = 1 Torr it can be well approximated by the following expression k 25 The pressure dependence is weak, and the rate constant may be considered pressure-independent in the temperature range where the contribution of this channel becomes significant. Much more effort is needed to tackle the (R24) channel.
Experimental values at pressure about 1 Torr and at room temperature In the pressure range of interest (p < 1.5 Torr) the rate constant k 24 is in the fall-off regime, therefore, it depends on p and on the collision partner (buffer gas). Some experimental determinations at 1 Torr in helium buffer gas include (6. [41]. Later, the smaller values received criticism in [42]. In the present study, we regard values around 1.1-1.3 × 10 −10 as the reference magnitude at room temperature. Data with different buffer gases can also be found in [42]. According to figures 3 and 4 of that study, the rate constant is about (1.7 ± 0.2) × 10 −10 in argon and 2.8 × 10 −10 in SF 6 . This trend is consistent with the concept that a heavier collision partner is more efficient in collisional stabilization of the intermediate complex. Furthermore, in [39] no significant difference has been observed at 5 Torr in helium buffer as compared with argon. In [43] no influence could be detected with argon addition to the SiH 4 /H 2 mixture even at a much lower pressure of 9 Pa. On the other hand, in [22] somewhat lower values in helium than in argon (by a factor of ∼1.5) have been determined below 1 Torr. This was explained by poor thermalization of SiH 2 radical in collisions with much lighter atoms.
T and p dependence Channel (R24) features both pressure and temperature dependencies. The T dependence is quite strong: the value at 500 K is smaller by a factor of 6.7 than that at 300 K [23]. The p dependence is also substantial. The rate constant varies by a factor of ∼3.0 between 0.1 and 1 Torr at 300 K and even stronger at 500 K. In the following, we explain how these values were derived. The rate constant at various conditions was determined in a comprehensive experimental study [42], and its magnitude demonstrates an impressive agreement with the results of a later theoretical calculation [23] over the temperature and pressure range.
While the results presented in [23,42] are for the pressure above 1 Torr, measurements in the pressure range p < 1 Torr at room temperature have been carried out in [22]. The determined rate constants and their pressure trend, given in figure 4 in that work, show a very good agreement with some earlier experimental results and with theoretical curve of [39]. This curve is reproduced in figure 1 as a blue line and is employed in our model as the reference dataset k 24 (p,300 K). It amounts to 1.25 × 10 −10 cm 3 /s at the reference pressure 1 Torr.
To our knowledge, the only experimental determination of the rate constant at non-room temperature, besides the studies by Becerra [42], was reported in [43]. The measurement yielded (2.3 ± 0.5) × 10 −10 at 500 K and 9 Pa. This value, however, is rather close to the data at room temperature, cf figure 4 in [22].

Preparing expressions for Quantemol input
To handle the T dependence, the data from [23] (see figure 10(a) therein) was employed. The dataset k 24 (p,300 K) was manipulated in order to generate smooth continuations to the curves at 500 K and 700 K. Since the pressure range is downside limited by 1 Torr, the data needs to be continued into the lower pressures. The pressure fits in [7] are based of the absolute values from [38,40] which most likely underestimate the rate constant, as stated above. Multiplying by a factor of about 2.5, the absolute values can be adjusted to get the desired magnitude at 1 Torr. Even so, it seems problematic to use the resulting expressions at pressures below 1 Torr. The reference curve (blue line in figure 1) has a more convex shape in this pressure range. As a result, the rate constant is massively underestimated when using the fits by Perrin.
In order to generate input data for simulations with Quantemol software, the obtained datasets were fitted by expressions linear in p, in which case the reaction can be represented as a sum of binary and ternary processes. Since the datasets are strongly curved as functions of p, for each of the three temperatures two fit lines were produced: for a higher pressure and for a lower pressure range. The resulting coefficients are interpolated as functions of temperature, so that the rate constant takes the form k(p,T) = a 1 (T) + a 2 (T) × p. This is shown in table 4.

Surface reactivities
To finalize the discussion on SiH x chemistry, a few remarks are needed on the surface reactivities found in table 1. The values of ß and s for SiH 3 are taken from [4]. The magnitude of the sticking coefficient for SiH 2 originates from an atomic-scale simulation study [10]. Therein, it was shown that the overall sticking coefficient for different H concentrations was approximately 70%. This value is close to the earlier estimates of 0.8 and also to an experimental determination 0.6 [43]. The sticking coefficient of the ions are taken as s(SiH x + ) = 1, following the surface chemistry model presented in [6]. The ion-induced secondary emission coefficient was set to 0.033, as determined in [9].
Outline of Si y H x (y = 2, 3) plasma chemistry In the literature, Si 2 H x radicals are less studied than SiH x . What makes it more complicated is the fact that each Si 2 H x radical may exist in the form of different isomers with varied reactivities. Regarding Si 2 H 4 radical, this is discussed in [14]. There is no way to do modelling without reasonable simplifying assumptions. Here we briefly describe our approach. The reasoning builds largely on results from [12]. Therein, three Si 2 H x radicals have been detected: Si 2 H 2 , Si 2 H 4 and Si 2 H 5 . The fact that they survived in significant amounts on their way to the surface suggests that they might be represented by low-reactive isomers. Thus, our model assumes that all the three radicals exist in the form of low-reactive isomers. Specifically, this means that they do not react with the background gases SiH 4 /H 2 and have moderate surface loss probability. Using the example of Si 2 H 4 as discussed in [14], a highly-reactive form of this radical (silylsilylene) is produced in chemical reactions but then undergo a fast unimolecular isomerization towards a low-reactive form (disilene). In that study, the sticking coefficient for this isomer was assumed 0.4. In our model, we assume the surface loss probability ß(Si 2 H 4 ) = 0.4, and the sticking s(Si 2 H 4 ) = 0.1. For Si 2 H 5 , the surface reactivities are assumed the same as for SiH 3 . For Si 2 H 2 , we use ß(Si 2 H 2 ) = 0.5, as suggested in [12], and s(Si 2 H 5 ) = 0.1 (0.3 was used in [12], but in their discussion the authors lean towards a lower value). Si 2 H 3 radical was also included in our model, although in a highly reactive form only. It must be produced in R27 by substantial amounts, because SiH is a product of the primary dissociation. The data in [22] suggests that the rate constants of R27 and R24 are similar. On the other hand, high reactivity of Si 2 H 3 results in enhanced consumption through R28 in the gas phase and high reactivity at the surface. Since Si 3 H x species are of minor relevance for the plasma-chemistry under present conditions, they are represented by a single lumped species. It effectively serves to close the system of rate equations for Si 1 -and Si 2 -chemistry. Within this lumped species, stable Si 3 H 8 molecule has the most dominant contribution, with the radicals forming an insignificant part.
Computing with Quantemol software Fluid simulations were performed using Quantemol Virtual Tool software [44]. This software package provides a GUI to the HPEM code developed by M Kushner [45]. It is capable of modelling 2D geometries with an option to choose either Cartesian or cylindrical coordinates. The mesh is defined in a way that all the numeric cells are the same and remain unchanged while simulation. The Boltzmann equation is solved for the electron energy Table 4. Coefficients for the approximation k(p,T) [cm 3 /s] = a 1 (T) + a 2 (T) * p[Torr], 300 K < T < 700 K. Within the indicated pressure ranges, the deviation remains less than 10% with respect to the reference curve k 24 (p, 300 K) and to the datasets obtained for 500 K and 700 K by the procedure described in the text. Setup for the radio frequency discharge in [12] Simulations for case (1) and case (2) are based on a mesh of a similar shape but different dimensioning. The 2D simulation domain was built following the brief outline of the geometric configuration found in [12]. It is characterized by a cylindrical symmetry with electrode radius 4.1 cm (diameter 8.2 cm) and substrate-toelectrode gap of 2 cm. The radical fluxes were calculated at the substrate center, corresponding to the centered position of the sampling orifice.
The radical densities at the substrate are calculated from the fluxes which were measured by TIMS (threshold ionization mass spectrometry) technique. The neutral radicals that pass through a small sampling orifice get ionized by electron impact. The use of near-threshold energy ensures that only the direct ionization occurs, i.e., each Si y H x + ion (y = 0, 1, 2) originates from the corresponding Si y H x radical. For the purposes of mass-resolved detection, the ions are further guided by the ion optics to the quadrupole mass spectrometer.

Results and discussion
Model validation at room temperature, cases (1) and (2) For the purpose of testing the correctness of a plasma chemistry model, it is instructive to compare its predictions with experimentally measured densities of chemical species. Determination of many radicals simultaneously is rarely found in literature [8,12]. For benchmarking, we have chosen work [12] which contains a comprehensive study of the whole set of SiH x and Si 2 H x radicals at different discharge conditions. The experiments were carried out at room temperature at varying dilution ratio and pressure. We used these parameters as the input to our model and ran simulations with the respective geometric configuration.
In figure 3 and figure 4 we compare the radical densities of our simulations with those reported in table 1 of [12]. In that study, the densities were calculated from the measured surface fluxes. In our simulations, Quantemol software computes the fluxes onto the wafer as part of the output. Similarly, we use these fluxes to infer the densities. This is done in the following manner. The Quantemol algorithm considers two components of a surface flux, the first one corresponding to the fluid average velocity, the second one arising from the thermal velocity of molecules. In our case, the latter component is largely predominant. Since we are interested in number densities, a ratio of fluxes times the inversed ratio of thermal velocities would yield the ratio of densities. Fluxes of all the radicals have been normalized to SiH 3 flux, making SiH 3 always appear as 100% on the diagram. Finally, they have been multiplied by the inversed ratio of thermal velocities, the velocity varying as / T m 2 .At a common temperature T for all the radicals, the relative density of, say, atomic hydrogen is by a factor of / m m SiH H 3 »5.6 smaller than its relative flux.
The first column with results in table 1 of [12] was acquired in pure silane as feed gas (dilution ratio R 0 = H 2 : SiH 4 = 0), pressure of 0.3 Torr and room temperature. This is what we call case (1). These results are compared with our simulations in figure 3, the silane depletion was about 10% in both cases.
All in all, a good agreement can be stated. The ordering of the most abundant radicals is correctly reproduced and none of them is missing. While Si-containing radicals contribute to deposition, the role of H atom is the opposite, i.e., etching. Thus, for a correct prediction of the deposition rate, it is necessary to accurately account for H flux. Getting back to the discussion on the branching ratio of dissociation reactions R1/R2, it is noted that the value used in our model R1/R2 = 0.46/0.26 produces a relative H content that agrees very good to the findings of [12]. Whereas R1/R2 = 0.17/0.83 generates a large amount of excess H that cannot be dealt away, e.g., by assuming a higher surface loss probability ß(H) which is already 1.0.   lower temperature and (b) high dilution ratio [13]. Thus, from figure 5 in [13] we can deduce that ß(H) ≈ 0.25 for R 0 = 40. The value used in our simulation was ß(H)=0.4, which gives better agreement than 0.3 or 0.5. Since the data represents the density near the surface, it is more sensitive to the loss probability than the average bulk value.
An interesting and novel result of the work by Horvath et al was eliciting a significant role of Si 2 H 2 among the other Si 2 H x species. This radical is rarely included into plasma chemistry models, though. Meanwhile, the experimental data does not show Si 2 H 4 radical. Most likely, its density lies slightly below the sensitivity threshold of 6%.
Model validation at T = 250°C, case (3) After validating the model at room temperature, we applied it to the MAiA capacitively coupled plasma source. Its model configuration is shown in figure 2. The x-dimension of 2D simulation domain was shorter than that of the real setup in order to save the computation time. We present the results for a typical set of parameters for PECVD depositions. The process temperature is around 250°C, which is much higher than in cases (1) and (2) discussed above. A silane dilution of R 0 = 4.5 at total pressure p=0.225 Torr is applied. This is case (3) of our validation.
In contrast to the diagrams above, figure 5 shows the bulk average values, not the values at the surface in the center. A striking feature of these results is a high amount of hydrogen H. It is even larger than in case with R 0 = 40 (figure 4), which is counterintuitive. Two effects are responsible for this. Firstly, due to the difference in the surface loss probabilities of H and SiH 3 , the proportion of H atoms would be lower if the densities at the surface were shown in figure 5 instead. Secondly, the surface losses of the atomic hydrogen are accompanied by enhanced volume losses in case of figure 4. This can be explained in the following way. H radical is produced by the electron-impact dissociation of SiH 4 as well as H 2 molecules. The losses are (a) due to the recombination on the surface and (b) due to H-abstraction reactions in the volume, most notably R20. The relative importance of volume/surface losses is characterized by the second Damköhler number (Da II ). Comparing the cases of figure 4 and figure 5, Da II is estimated to be about 20 times larger in the first case, considering the differences in temperature, geometry, total and partial pressures and surface loss probability. This provides an indicator for the increased efficiency of the volumetric losses.
The relative number density of SiH 2 is higher than on the previous plots due to the temperature dependence of reaction R5, which leads to reduced losses of this radical at a higher temperature.
The resulting spatial distributions of plasma species typically show maximum in the region under the edge of the powered electrode. This is illustrated by a 2D mapping of SiH 3 radical in figure 6.
The indicated power density is reduced to the dimension of the powered electrode. In the simulation, the pumping residence time of a molecule was 0.3 s. This is close to the experimental value as estimated from the chamber dimensions, gas inflow and pumping. The residence time and the dissipated power determine the resulting gas depletion. With plasma on, SiH 4 density was only 24% of the initial plasma-off condition. Due to this high degree of depletion, the influent dilution ratio R 0 = 4.5 of the feed gases turns into the actual value R = 22.2 in the operated plasma chamber.
Silane is known to produce electronegative plasmas, i.e., large amounts of negative ions. In our case, the simulated number density of SiH 3 − is by a factor of 1.5 larger than electrons. The literature on experimental determination of electronegativity in silane-containing plasmas is scarce. Thus, in [47] the number density of the negative ions was estimated about 19 times larger than that of electrons for the case of R 0 = 1 dilution. In [48], with silane content of 0.5% (diluted in He) the densities of the negative ions and the electrons were found to be of comparable magnitude. When modelling a real experimental situation, it has to be made sure that the dissipated power matches. Indeed, in both the experiment and the simulation the dissipated power density was around 0.15 W cm −3 . This value is reduced to the width of the powered electrode. An efficient path for power dissipation is the gas heating due to electron-impact rovibrational excitation. These reactions are not explicitly present in the plasma chemistry set. However, in Quantemol software, the power losses due to excitation of the four vibrational modes of SiH 4 molecule are automatically taken into account when solving Boltzmann equation. The problem of power dissipation in H 2 /SiH 4 mixtures has also been tackled in [5]. The higher the H 2 content, the less efficient this dissipation channel is. From figure 3 in [5] we learn that, with actual R = 22.2, the power dissipated through this channel is about 12%.
The thickness of the deposited layer at this condition was measured by means of ellipsometry. It is about 18 nm for the inline process. Considering the speed of the tray 55 cm s −1 and the width of the plasma source about 12 cm, this thickness translates into deposition rate of ≈ 83 nm min −1 .
This value can be compared to the one obtained via simulation. The deposition rate profile is shown in figure 7. It was computed using the surface fluxes of radicals and ions and their sticking coefficients given in table 1. Also included is etching by H atoms, using the temperature-dependent expression for etch probability found in [8]. The non-uniform spatial profile of the deposition rate was averaged over the wafer, in the way like it is smoothed out by the inline motion of the carrier in the real setup. The resulting value of 40,6 nm min −1 was then reduced to the width of the powered electrode. This finally yields the deposition rate of 87 nm min −1 , which is in a very good agreement with the one from the experiment.