COMSOL 6.3 - Surface Chemistry Tutorial Using the Plasma Module

Surface Chemistry Tutorial Using the Plasma Module

Surface chemistry is often the most important and most overlooked aspect of reacting flow modeling. Surface rate expressions can be hard to find or not even exist at all. Often it is preferable to use sticking coefficients to describe surface reactions because they can be estimated intuitively.

Model Definition

The tutorial model simulates outgassing from a wafer during a chemical vapor deposition (CVD) process. Careful attention is paid to the overall mass balance in the system and the difference between the mass averaged velocity and diffusion velocity is explored.

The same physical problem is investigated first with a global (volume-averaged) model and after with a space dependent model. An initial study of a problem with a global model is a good approach when there is limited knowledge about the system. One reason is that the global model has reduced computational times, allowing for fast sweeps over parametric spaces. The fast computational time of global models is also a great advantage to investigate models with complex chemistries.

In some cases a volume-averaged model might even be good enough to extract the desired information from the model. If, however, space dependent information is necessary, a space dependent model can be implemented in a subsequent study.

The geometry and operating principle for the model is shown in Figure 1. Initially the closed container is full of 99.9% silicon hydride and 0.1% hydrogen (measured by molar content).

Figure 1: Geometry and basic operating principle for the surface chemistry tutorial model.

A surface reaction begins to occur on the wafer which consumes the silicon hydride, releases hydrogen into the domain, alters the composition of the absorbed species on the wafer surface and deposits bulk silicon. The following reactions are considered on the surface of the wafer:

Table 1: surface reactions considered.
reaction	Sticking coefficient or rate constant
SiH4+2Si(s)=>Si(b)+2SiH(s)+H2	2E-4
SiH(s)=>Si(s)+0.5H2	3 1/s

The (s) here denotes “surface species” which means that the species only exists on the surfaces where the reaction is occurring. To indicate that a species is bulk, append (b) to the end of the species. Since bulk species cannot participate in surface reactions, they must only be products and not reactants in a surface reaction. The net result of these two competing reactions is SiH4=>Si(b)+2H2, which means that it is expected that the silane is replaced by hydrogen inside the reactor, and layers of silicon deposited on the wafer surface.

model equations — surface reactions and surface species

The surface reaction rate for reaction i is given by

where ck is the molar concentration of species k. The rate constant kf can be given by

where Γtot is the total surface site concentration (SI unit: mol/m2), m is the reaction order minus 1, T is the surface temperature, R is the gas constant, Mk is the molecular weight, and γi is the sticking coefficient. When defining a surface reaction that contains no gas-phase reactants consider specifying the reaction using a rate coefficient. Otherwise, if using the equation above the mean gas-phase molecular weight is used.

For the surface species the following equations are solved:

where Zk is the site fraction (dimensionless), Rsurf,k is the surface rate expression (SI unit: mol/m2). The quantity Zk is the surface equivalent of the mole fraction on the volumetric level. That is, the sum of the site fractions must always equal one:

Empty sites are not accounted for in reaction rate coefficients, so all species which can occupy a surface must be included in the model. The surface site concentration is given by

For each of the bulk surface species, the following equation is solved for the deposition height:

where hk is the total growth height (SI unit: m), Mk is the molecular weight (SI unit: kg/mol), and ρ is the density of the bulk species (SI unit: kg/m3).

domain equations

Inside the domain, the Navier–Stokes equations are solved for the fluid velocity. The mass fraction of hydrogen is computed by solving:

where w is the mass fraction of hydrogen. For detailed information on the transport of the nonelectron species see Theory for the Heavy Species Transport Interface in the Plasma Module User’s Guide. The mass fraction of silane is not directly computed. Its value comes from the fact that the sum of the mass fractions must equal one. The gas temperature is computed by solving the energy equation.

boundary conditions

Fluid Flow

The following boundary conditions are used. The mass averaged velocity is constrained using:

where Mf is the inward (or outward in this example) mass flux which is defined, from the surface chemistry as:

where sk is the surface rate expression for each species which comes from summing the surface reaction rates multiplied by their stoichiometric coefficients over all surface reactions:

Hydrogen Mass Fraction

The flux of hydrogen at the surface of the wafer comes from the surface reactions:

Energy Equation

For the energy equation, the following boundary condition is used:

where hi is the molar enthalpy change due to reaction i.

Global model Equations

The global model used in this work considers that the spatial information of the different quantities in the plasma reactor can be treated as uniform. Without the spatial derivatives the numerical solution of the equation set becomes considerably simpler and the computational time is reduced. For detailed information on the global model see Theory for the Heavy Species Transport Interface in the Plasma Module User’s Guide. For heavy species, the following equation is solved for the mass fraction:

where ρ is the mass density (SI unit: kg/m3), wk is the mass fraction, and Rk is the rate expression (SI unit: kg/(m3.s)). The fourth term on the right hand side accounts for surface losses and creation, where Al is the surface area, hl is a dimensionless correction, V is the reactor volume, Mk is the species molar mass (SI unit: kg/mol) and Rsurf,k,l is the surface rate expression (SI unit: mol/(m2.s)) at a surface l. The last term on the right-hand side is introduced because the species mass balance equations are written in the nonconservative form and it used the mass-continuity equation to replace for the mass density time derivative. In the last term, Mf,l is the inward mass flux of surface l (SI unit: kg/(m2·s)). The sum in the last two terms is over all surfaces where there are surface reactions.

To take into account possible variations of the system total mass or pressure, the mass-continuity equation is also solved:

The gas temperature is set to a constant value. The model for the surface chemistry is exactly the same.

Results and Discussion

In this section, both space-dependent and global model results are presented and discussed. For this particular physical system, the results from the global model and the space dependent model are in a very good agreement.

The y-component of the mass-averaged velocity is plotted in Figure 2. The mass-averaged velocity is negative at the surface of the wafer. This means that overall, mass is leaving the system. This is to be expected since silane is being consumed, which has a molecular weight of 0.032 kg/mol and replacing it with hydrogen, which has a lower molecular weight of 0.002 kg/mol.

Figure 2: Plot of the mass averaged velocity field after 5 seconds.

The surface reactions begin to consume the silane, which results in concentration gradients within the reactor. The outward mass flux at the wafer surface leads to a mass-averaged velocity everywhere inside the reactor. The combination of concentration gradients and convection due to the mass-averaged velocity tends to draw silane toward the wafer. The first surface reaction is exothermic while the second surface reaction is endothermic. The amount of heat released by the exothermic reaction dominates, so the temperature begins to increase. The temperature is plotted in Figure 3 and is highest on the wafer surface.

Figure 3: Plot of the change in gas temperature after 5 seconds. The higher temperature is observed at the wafer surface due to the heat released in the surface reactions.

From the net reaction SiH4=>Si(b)+2H2 there should be a molar inflow of hydrogen at twice the rate that silane is leaving. Although this condition is not applied explicitly, it is implicit from the equations solved in the Heavy Species transport interface. Despite the fact that there is a negative mass-averaged velocity, there is a positive diffusion velocity for the hydrogen at the wafer surface. The y-component of the diffusion velocity for hydrogen is plotted in Figure 4.

Figure 4: Plot of the y-component of the diffusion velocity for hydrogen after 5 seconds. Even though the mass-averaged velocity is negative, the diffusion velocity is positive, indicating an inflow of hydrogen and a net outflow of total mass.

Figure 5 plots the integrated pressure in the reactor divided by the initial average pressure in the reactor. Once the problem reaches steady state, the average pressure in the reactor has increased by a factor of 2. This is expected since every mole of silane is being replaced by two moles of hydrogen. The silane in the reactor is completely consumed after around 200 seconds.

Figure 5: Plot of the ratio of the average pressure by the average initial pressure computed with the global model (blue) and the space dependent model (green).

In Figure 6, the total mass initially present in the reactor divided by the total mass in the reactor is plotted. The total mass in the reactor drops by a factor of 8, which is expected since the silane with molecular mass 0.032 kg/mol is replaced with two moles of hydrogen, which has a molecular weight of 2·0.002 = 0.004 kg/mol.

Figure 6: Plot of the ratio of the total mass by the total mass computed with the global model (blue) and the space dependent model (green).

Another way of verifying the correctness of the model is to compare the mass lost inside the reactor to the mass accumulated on the wafer surface. The two quantities should be equal. The total gain in mass from surface and bulk species is compared to the total loss in mass from the reactor in Figure 7. The two curves agree very well, indicating that the total mass in the entire system is conserved.

Figure 7: Plot of the total mass lost in the domain (blue: global model, and red: space dependent model) and the total mass gained due to Silicon deposition on the wafer surface (green: global model and cyan: space dependent model).

The ultimate goal of most CVD models is to determine the total growth height and growth rate on the surface of the wafer. The total growth height is plotted in Figure 8 and saturates at about 159 angstroms. The total growth rate saturates because all the silane in the reactor is consumed after around 200 seconds. The accumulated growth height is also very uniform across the surface of the wafer. The global model and space dependent model results for the total growth height agrees very well (not shown).