Analysis of hydrogen plasma in a microwave plasma chemical vapor deposition reactor

Carbon nanostructures such as carbon nanotubes (CNTs), graphitic sheets and nano-petals, and graphene find a wide variety of applications in the fields of engineering and bio-sciences owing to their exceptional mechanical,^1,2 thermal,³ and electrical⁴ properties. Due to the rapid increase in their utilization over the last decade, there is a need for synthesizing high quality carbon nanostructures at reduced cost. Chemical vapor deposition (CVD) is a method for producing such nanostructures on a large scale by depositing gaseous reactants onto a substrate that is placed in a reactor at ambient temperature. The precursor gas molecules dissociate in the reactor and recombine to create a material film on the substrate. The substrate temperature is a key parameter for the growth of the material. The products of CVD are of high quality, hardness, purity, and density⁵ when compared to materials produced by other processes.

Microwave plasma chemical vapor deposition (MPCVD) is one of the commonly used CVD processes. The first step in this process is the introduction of pure hydrogen into the reactor. The electric field generated by microwaves ionizes the hydrogen, igniting plasma in the reactor. The next step is the introduction of methane after the reactor reaches a steady state. Methane dissociates in the presence of plasma resulting in the deposition of carbon on the substrate for the growth of various carbon nanostructures. A schematic of a methane based plasma assisted CVD process is as shown in Fig. 1. The MPCVD technique is energy efficient, allows for low temperature synthesis⁶ of carbon nanostructures, and produces stable plasma that can be sustained over a long period of time. MPCVD systems used for diamond deposition have been modeled previously to study pure hydrogen discharges^7–14 as well as hydrocarbon discharges.^7,15,16 Fluid models have also been used to investigate radio frequency discharges used for plasma assisted CVD processes.^17,18

The growth rate and quality of the carbon nanostructures depend on the number density of atomic hydrogen and the concentration of certain precursors on the surface of the substrate, which in turn depend on the discharge and operating conditions.¹⁹ H₂ is the major species in the H₂/CH₄ plasma for the MPCVD process. Although the introduction of small amounts of hydrocarbons in the feed gas may affect the plasma interaction with the substrate, the main features of the plasma discharge configuration can be obtained by studying pure hydrogen plasma inside the reactor.⁷ In the present work, we model such pure hydrogen plasma and its interaction with the gas, because information on the electric field, electron number density, electron temperature, and gas temperature are indispensable for understanding and optimizing the operation of MPCVD systems. The electric field is modeled by a high frequency Maxwell solver. This electric field is used to determine the electron number density by Füner's model^10,11,20 and electron temperature based on the model by Yamada et al.¹¹ The gas temperature is calculated by a heat transfer solver. The various solvers and models are coupled using COMSOL Multiphysics. A similar model has previously been demonstrated by Huang.²¹ A Boltzmann equation solver for electrons in weakly ionized gases under given local electric field, called BOLSIG+,²² is used in conjunction with these solvers. The aim of this work is to determine hydrogen plasma characteristics using simple modeling techniques, and the aforementioned solvers and models are chosen because they can be easily coupled without detailed chemistry modeling. The modeled gas temperatures are compared to experimentally measured temperatures of hydrogen using optical emission spectroscopy (OES) techniques.

The MPCVD reactor, a Seki Technotron Corp. AX5200S model, has been used extensively to grow carbon nanostructures such as CNTs,²³ graphene,²⁴ and contiguous graphitic fins²⁵ over a variety of substrates. The system is shown in Fig. 2. Microwaves are generated by a 1.5 kW (2.45 GHz) ASTeX AX2100 microwave generator and a 3.5 kW radio frequency power supply. The microwave field passes through a rectangular waveguide, a TE-TM mode converter and into the MPCVD reactor which is isolated from the rest of the system by a quartz plate. The reactor or plasma chamber includes a graphitic susceptor that can move along a vertical axis and a molybdenum puck which is used to concentrate the plasma near the substrate. In the case of graphitic nanopetal growth, a dielectric pillar is placed on the puck to concentrate the plasma.

The rotational temperatures of the Q-branch of the H₂ $d3Πu→a3Σu+$ (0,0) Fulcher (α) band are measured using OES. The experimental setup is shown in Fig. 3, and an example of the measured H₂ spectrum is shown in Fig. 4. The emission from the plasma is collected and transmitted using an optic fiber cable to a spectrometer (Acton SP-2756 Princeton Instruments). A CCD camera (PIXIS 256E Princeton Instruments) is used to detect the light. The spectrometer is calibrated using dual mercury and neon-argon lamps in the wavelength range from 300 to 1000 nm.

The intensity (I) of the emission from an upper rotational level (⁠$EJ′$⁠) to a lower rotational level (⁠$EJ″$⁠) is given by

where N_u is the population of the excited state, A_ul is the transition probability (or Einstein coefficient), $ν̃$ is the wavenumber, and h and c are the Planck's constant and speed of light, respectively. For a rotational transition with a particular electronic state and vibrational band that follows Boltzmann distribution, the emission intensity can be rewritten as

where $SJ′J″$ is the Hönl-London factor which is taken from Qing et al.,²⁶ λ is the wavelength of the transition obtained from Crosswhite and Dieke,²⁷ k_B is the Boltzmann constant, $EJ′$ is the rotational energy of the upper excited rotational state, and T_rot is the rotational temperature of the upper state (d³Π_u). The rotational energy (⁠$EJ′$⁠) is calculated assuming rigid rotor model approximation

where $J′$ is the rotational quantum number for the upper energy level and B_v is the rotational constant obtained from Qing et al.²⁶ Table I shows the Hönl-London factors, rotational energies, and wavelengths²⁷ for Q₁, Q₂, Q₄, Q₆, and Q₈. These wavelengths are selected in the Boltzmann plot since they have no overlap with other transitions.²⁸ 

An example of the measured H₂ spectrum with indication to the selected rotational lines is shown in Fig. 4. The Boltzmann plot is presented as an inset for this spectrum where $ln(Iλ4/SJ′J″)$ is plotted as a function of the corresponding rotational energy (⁠$EJ′$⁠). The rotational temperature can be deduced from the slope of the plot (⁠$−1/kBTrot$⁠). In the current study, the rotational temperature of the upper state of H₂ (d³Π_u) is assumed to be the same as the H₂ translational state. Comparison with previous temperature measurements using CARS in the same system supports the above assumption within experimental uncertainty of 10%.²⁹ 

A Princeton PI-MAX 4 iCCD camera is used to obtain optical emission images of the plasma for qualitative comparison with the modeling results.

The first step in modeling the plasma in the MPCVD reactor is to determine the electric field because the electron number density and electron temperature depend directly on it. The electromagnetic (EM) simulations treat Maxwell's equations in the form

where $k02=ω2ϵ0μ0$⁠, ω is the angular frequency of the microwave field, $E→$ is the electric field vector which consists of the amplitude as well as the time component, ϵ₀ and ϵ_r are the vacuum and relative permittivities, μ₀ and μ_r are the vacuum and relative permeabilities, and σ is the conductivity. For the first step, vacuum conditions are assumed and a solution is obtained for the electromagnetic field. However, the presence of plasma in the domain changes the material properties of the system by increasing the conductivity of the gas. The plasma properties are evaluated as^11,20

where ω_p is the plasma frequency given by

Here, e represents the charge of an electron, m is the mass of an electron, n_e is the number density of electrons, and ν_m is the collision frequency between electrons and other species evaluated as

where N is the number density of the gas, v is the mean electron velocity, and σ_s is the collisional cross section. For pure hydrogen plasmas, this collision frequency can be approximated as^10,11,20

where p_g is the gas pressure in Pa and T_g is the gas temperature in K. These plasma properties are updated in the following iterations by the EM solver.

The electron number density is evaluated based on Füner's model^9,10,14 as

where E_m is the minimum electric field strength needed to sustain a discharge, which is taken to be half of the breakdown voltage.³⁰ From MacDonald et al.,^31,32 we find that its value is

Because the plasma is generated by microwaves, the radiation must be able to penetrate into the plasma in order to ionize molecules. This happens only when $ωp2≤ω2+νm2$⁠. However, the skin depth for the plasma is on the order of the size of the plasma region. Therefore, the maximum electron number density constraint is not imposed, but is used to determine the minimum electron number density. The value of n_e,min that we use in the model must be lower than the value of n_e,max. The maximum electron number density in the plasma is thus given by

where ν_m is given by Eq. (9). At the two pressures, we have

Based on this constraint as well as experimental electron number density measurements on similar systems,^12–14,33 we use $ne,min=1×1017m−3$⁠. The value of γ is calibrated to match the gas temperature predicted by the model with experimental measurements at 500 W and 30 Torr. The variation of predicted gas temperatures with γ is shown in Table II. Since the value at $γ=3×1013V−1m−2$ matches best with the experimental measurement of 1441 K, we choose this value of γ for all computations.

We use the model given by Yamada et al.¹¹ to determine the electron temperature. Simplifying the conservation of energy for electrons generated by microwaves, the rate of energy transferred from the microwaves to the electrons per unit volume, Q_MW, is shown to equal the rate of energy loss of electrons during collisions with heavy species, Q_c. The collisional power loss term is given by

where E_s, n_s, and k_s are the threshold energy, species number density, and reaction rate coefficient for reaction s. The term Q_c/n_e is approximated by a function given by

where E_L, p_L, and α_L are obtained by fitting a curve over a range of electron temperatures to match the two sides of the following equation:

For the present study, the reactions considered are as shown in Table III.

The rate coefficients for excitation and ionization are obtained from the Boltzmann solver²² with the cross section data taken from the LXCat database.^34–36 The rate coefficient for dissociation is obtained from the Arrhenius equation given by³⁷ 

where A_d is the pre-exponential factor equal to $1×10−14m3s−1$⁠. Hydrogen undergoes a number of excitation reactions. The variation of Q_c/n_e with electron temperature for the excitation reactions with the highest values of Q_c/n_e in the range of electron temperatures that are of importance to the current simulations is as shown in Fig. 5. We choose the excitation reaction with E_s of 8.9 eV because at the maximum of electron temperatures seen in the reactor, this reaction has the highest contribution to the collision term. The reaction with E_s of 0.52 eV has a greater contribution at low electron temperatures. The error in the estimation of Q_c/n_e due to the omission of this reaction at T_e = 2 eV is 60% at 10 Torr and 25% at 30 Torr.

The curve fitting parameters used to determine the electron temperature are given in Table IV. The approximation, (Q_c/n_e)_app, agrees well with the actual collision term, Q_c/n_e, as indicated by Fig. 6.

By approximating the collision term and setting it equal to the microwave power density absorbed, the expression obtained for electron temperature is given by¹¹ 

where

A heat transfer solver is used to determine the gas temperature, T_g, using the 2D heat conduction equation with a volumetric heat source term, Q_g given by

where k is the thermal conductivity of the gas. The term Q_g represents the power density imparted to the neutral gas molecules when they collide with the free electrons in the plasma chamber. Here, we consider only hydrogen molecules. The chamber is a cylindrical shell. For a gas temperature of 2500 K and ambient temperature of 300 K, the Rayleigh number for a vertical cylinder of height 130 mm is found to be 250, which is much lower than the critical value of 3800 for a similar configuration.³⁸ Therefore, thermal convection is not taken into account in our model. The heat source term³⁹ is given by

where δ is a fractional energy loss term which is a measure of the fraction of energy transferred from electrons to neutral gas during collisions and k_B is the Boltzmann constant. Using the approximate electric field obtained from the EM solver and the gas temperature, the total power loss is calculated by the Boltzmann solver²² for various reactions of hydrogen. The term δ is then found using the following equation:

where dϵ_e/dt is the total power loss which is the sum of elastic and inelastic power losses and ϵ_e is the mean energy of electrons. The variation of δ with the reduced electric field for rotational, vibrational, and electronic excitation reactions of H₂, ionization of H₂, and elastic collisions between electrons and H₂ molecules at a temperature of 1000 K is shown in Fig. 7. The δ value with all the processes included is used to calculate the gas temperature.

The various solvers are coupled as shown in Fig. 8. The Maxwell solver determines the electric field, which is used by the user defined functions (UDFs) along with the properties from the Boltzmann solver to determine the electron number density, plasma material properties, and electron temperature. They are used to calculate the heat source term for the heat transfer solver, which determines the gas temperature. Iterations are carried out till convergence of residuals is obtained.

The plasma region can be modeled as a 2D axisymmetric region owing to the fact that it is cylindrical. A part of the transmission line, the converter, the quartz plate, and the plasma chamber are modeled axisymmetrically. The susceptor is modeled at the bottom of the plasma chamber, with the puck placed concentrically above its top surface as shown in Fig. 9. The grey region is the part modeled for the EM simulations and the green is modeled for both EM and heat transfer simulations. In some simulations, the susceptor and puck assembly is raised or lowered by a few mm in the reactor to study the effect of the susceptor position on the plasma characteristics. For other simulations, an axisymmetric pillar is modeled above the puck in order to study the effect of change in the geometry of the susceptor stage.

The boundary conditions applied are shown in Fig. 9. For the pure EM simulations without the plasma effects included, vacuum is assigned to the plasma reactor. When the EM solver is coupled to the plasma and heat transfer models, the reactor is composed of an inlet, outlet, walls, and susceptor-puck assembly. For the heat transfer solver, the temperature of the susceptor-puck assembly is maintained at experimentally observed values. A thermocouple placed below the susceptor surface is used to obtain these measurements while running the system with hydrogen plasma under various conditions. The initial temperature of gas in the plasma chamber is taken to be 293 K. The pillar is modeled as a thermal insulator.

A physics controlled mesh is used in COMSOL to discretize the domain with greater refinement in corner regions. A grid convergence study is conducted with different grid sizes and input conditions to ensure convergence of gas temperature within the plasma chamber. The results of the study for an input power of 500 W and gas pressure of 30 Torr are shown in Fig. 10. When the number of elements is increased beyond 9304, it does not change the temperature by a significant amount and so this grid has been adopted for the present calculations.

When microwave input power is increased at constant pressure, the electric field inside the plasma chamber increases because of an increase in the energy density of the system. An increase in power always leads to the expansion of the high field region as seen in Fig. 11(a). At 10 Torr, when the power is increased, the maximum electric field reached in the plasma chamber reduces, but the average field increases due to this expansion. At 500 W and 10 Torr, a residual high field region appears below the quartz plate which may lead to the formation of a two-ball plasma if the power is increased further. Such a configuration is not desirable because it not only takes away power from the region of interest but may also cause structural damage to the system. An increase in the pressure of the gas in the chamber leads to an increase in the collision rate, thereby increasing the energy loss inside the system. This reduces the number density of electrons and thus increases plasma permittivity which in turn is fed into the EM solver and results in a decrease in the electric field. Higher pressure also causes the high electric field to be concentrated in a smaller region as seen in Fig. 11(b). Higher values of maximum field are reached at higher pressure, but the average field in the chamber is lower because the high field region shrinks.

The electron number density depends directly on the electric field as given by Eq. (10). As the power is increased at constant pressure, the average electron number density increases, and as the pressure is increased, it decreases. The effects of input power and gas pressure on electron number density inside the plasma chamber are shown in Fig. 12. The increase in power goes into increasing the volume of the plasma region. At 10 Torr, the plasma region is slightly detached from the puck. Previous work reported by Hassouni et al.⁷ on similar systems demonstrated the change in electron number density while keeping the plasma volume a constant. This requires a simultaneous change in the power and pressure. However, in the present work, the power is kept constant as pressure is varied and increasing the pressure is seen to concentrate the plasma in a smaller region, thereby reducing the average electron number density in the chamber. The electron number density values stay fairly uniform both axially as well as radially within the plasma region.

The electron temperature depends directly on the electric field and inversely on the gas pressure. Therefore, it is seen to increase with an increase in the power and decrease with an increase in pressure which is seen in Fig. 13. The electron temperature stays fairly constant within the plasma region axially as well as radially.

The gas temperature depends on the energy transferred to gas molecules by electrons during collisions. The electron temperature increases as the input power is increased at constant pressure. This makes more energy available to be transferred thereby increasing the gas temperature. An increase in the gas pressure reduces the electron temperature, but increases the number of collisions. Hence, more energy is transferred and so gas temperature increases. The average gas temperature, however, reduces with increasing pressure on account of the plasma region being smaller. The effect of gas pressure and input power on gas temperature can be seen in Fig. 14.

The experimentally obtained temperatures have been compared to the modeling results for the case without the pillar. Fig. 15 shows the comparison of these temperatures at a point 5 mm vertically above the center of the puck. This point is indicated by a star in Fig. 9. The comparison indicates that the computational values agree well with the experimental observations, especially at higher input powers and gas pressures. At 10 Torr, the numerical results are higher than the experimental measurements by approximately 150 K. For both sets of results, the gas temperature has a weak dependence on power and a very strong dependence on pressure. On an average, the gas temperature at 5 mm above the center of the puck increases from about 1100 K at 10 Torr to about 1500 K at 30 Torr.

The contours of electric field at 500 W and 30 Torr are shown in Fig. 16(a). When the pillar is included, the electric field reaches very high values close to the top of the pillar. These values drop rapidly through the high field region both axially as well as radially. Consequently, even though maximum field values reached are very high with the pillar, the average field values in the chamber are very similar to when the pillar is not included. The contours of electron number density at 500 W and 30 Torr are shown in Fig. 16(b). The inclusion of the pillar concentrates the plasma region close to the top of the pillar. Very high values of n_e are reached close to the top of the pillar and the values fall rapidly through the plasma region both axially as well as radially. The contours of electron temperature at 500 W and 30 Torr are shown in Fig. 16(c). When the pillar is included, very high values of electron temperature are reached close to the top of the pillar, and the value decreases throughout the plasma in the axial and radial directions.

Fig. 16(d) shows the gas temperature in the plasma chamber at 500 W and 30 Torr. High electron number density close to the top of the pillar increases the number of collisions and high electron temperature in this region causes more energy to be transferred during each collision. Consequently, the gas temperature reaches high values close to the top of the pillar when it is included. This value decreases throughout the plasma region in the axial and radial directions as seen in Figs. 17 and 18. When the pillar is not included, the maximum value of gas temperature is reached at about 20 mm above the puck surface and reduces beyond the maximum point in the axial direction as shown in Fig. 17. Radially, the gas temperature stays fairly constant through the plasma as shown in Fig. 18.

Experimental images of plasma's optical emission in the range 200 to 800 nm (Princeton, PI-MAX 4 iCCD with 500 ms gate setting) have been captured at 500 W and 30 Torr with and without a pillar of 9.6 mm height. This pillar height is by necessity slightly less than that discussed above to allow full visual access through the available viewport. The resulting images shown in Fig. 19 are qualitatively similar to those of the model (cf. the predicted cross sections of various plasma-affected fields in Fig. 16). The plasma intensity is clearly higher and more concentrated around the pillar, which also elevates the plasma region vertically (axially) while contracting it horizontally (radially). Many species in the plasma, such as H₂ and H, exist in excited electronic states due to electron impact excitation^40,41 and consequently, the measured emission intensities corroborate the shapes of electron number densities predicted by the model with and without the pillar.

In our model for electron temperature, we assume that all the power absorbed from the microwaves by the electrons is transferred to the neutral gas during collisions. The power absorbed from the microwaves may be calculated as $QMW=e22mνm(ω2+νm2)ne|E→|2$⁠. The efficiency of power absorption is given as (Q_MW/input microwave power) × 100. Fig. 20 indicates that both these quantities are higher at low pressure and the inclusion of the pillar improves them at high pressure. This is because the pillar concentrates plasma of high electron number density and electron temperature close to the top, and the power is utilized in reaching higher values of electron number density and temperatures than in expanding the plasma.

The susceptor stage is moved without the inclusion of the pillar at 500 W and 30 Torr. The condition used thus far is used as the reference position (Y_sus = 0 mm). When the stage is lowered by 20 mm (Y_sus = −20 mm), the region of high electric field moves to the center of the plasma chamber by detaching from the puck. This causes the plasma to move away from the region of interest as can be seen in Fig. 21. Even though very high values of gas temperatures are reached in this case as shown in Fig. 22, they are not useful for the growth of nanostructures.

When the susceptor is raised by 25 mm (Y_sus = 25 mm), the plasma reaches higher values of electric field, electron number density, and electron temperature close to the puck. The gas temperature reduces by 3%. The plasma region stays attached to the puck. These conditions are favorable for deposition. However, a residual field appears to be forming near the top of the chamber which may lead to the formation of two-ball plasma if the input power is increased. When the susceptor is raised even more to Y_sus = 50 mm, we see the formation of secondary high field and plasma regions on the sides of the susceptor walls. This takes away power from the plasma and inhibits growth near the substrate. The secondary plasma may also damage the sides of the susceptor due to which this configuration is not desirable.

The numerical results with changes in susceptor position and inclusion of the pillar are not compared to experimental measurements quantitatively because the beam optics involved make it challenging to conduct such experiments.

A numerical model has been built to analyze hydrogen plasmas in an MPCVD system using the high frequency Maxwell solver and heat transfer solver in COMSOL Multiphysics. These solvers have been coupled with user defined functions as well as the Boltzmann equation solver, BOLSIG+, to determine plasma properties. Variation of discharge characteristics with change in input properties, susceptor stage geometry, and susceptor position has been studied. The gas temperatures obtained from modeling have been compared to experimental measurements of rotational temperatures of hydrogen in the same system, and the two show good agreement, especially at higher powers and pressures.

The simulation results indicate that an increase in the input microwave power leads to an increase in the electric field, electron number density, electron temperature, as well as gas temperature in the plasma chamber, and is accompanied by an expansion of the plasma. An increase in gas pressure reduces the electric field, electron number density, and electron temperature, and increases the gas temperature. The microwave power absorbed is higher at lower pressures. When the pillar is not included, the electric field, electron number density, and electron temperature are found to stay almost constant within the plasma region both axially and radially. The gas temperature almost stays constant radially, but in the axial direction, it peaks at about 20 mm above the puck surface. Both the numerical and experimental results show a weak dependence on input power and a strong dependence on gas pressure.

The inclusion of the pillar above the surface of the puck causes high electric field to be concentrated near its top and thereby, at 500 W and 30 Torr, the maximum electric field increases by 50%. This increases the maximum value of electron temperature by 70%, the maximum gas temperature by 50%, and the maximum electron number density by about 70%. These high values are seen right above the top of the pillar's topmost surface and these values decrease gradually through the plasma both in the axial as well as the radial directions. The electron number density contours show qualitative agreement with optical emission images of the plasma at similar operating conditions with and without the pillar. The inclusion of the pillar is favorable for growth of nanostructures because it concentrates the plasma of high strength in a smaller region near the substrate and improves the efficiency of microwave power absorption.

When the susceptor position is lowered, the plasma moves from the puck surface to the center of the plasma chamber, away from where the substrate would be placed which is not favorable to nanostructure growth. When the susceptor is raised by a small value (25 mm), it causes the plasma to be concentrated in a smaller region above the puck surface, increasing the values of electric field, electron number density, and electron temperature which is favorable. However, if the power were to increase further, the configuration may produce two ball plasma which is highly undesirable. When the susceptor is raised by a large value (50 mm), it causes concentration of plasma on the sides of the stage which causes a large part of the power to be taken away from the main plasma region above the puck surface, which is also not desirable. This indicates that the best position for the susceptor would be to raise it by less than 25 mm from its initial position.

This work was supported by NSF CMMI Grant No. 1344654. We thank Di Huang for his work on the initial modeling of the system and Dr. Abbas Semnani for his insight on the TE-TM mode conversion. We also thank Prof. Sergey Macheret for his suggestion of the maximum electron number density.