Model of PE-CVD Apparatus: Verification and Simulations

Abstract

In this paper, we present a simulation of chemical vapor deposition with metallic bipolar plates. In chemical vapor deposition, a delicate optimization between temperature, pressure and plasma power is important to obtain homogeneous deposition. The aim is to reduce the number of real-life experiments in a given CVD plasma reactor. Based on the large physical parameter space, there are a hugh number of possible experiments. A detailed study of the physical experiments in a CVD plasma reactor allows to reduce the problem to an approximate mathematical model, which is the underlying transport-reaction model. Significant regions of the CVD apparatus are approximated and physical parameters are transferred to the mathematical parameters. Such an approximation reduces the mathematical parameter space to a realistic number of numerical experiments. The numerical results are discussed with physical experiments to give a valid model for the assumed growth and we could reduce expensive physical experiments.

1. Introduction

We motivate our study by simulating the growth of a thin film by PE-CVD (plasma-enhanced chemical vapor deposition) processes; see [1, 2]. Such technical processes are very complex and real-life experiments are enormously extensive and expensive. Based on a large physical parameter space, the number of possible experiments involves at least the variation of all possible parameters. Such a large number of experiments can be reduced with the help of numerical experiments based on a mathematical model. Such modelling results are based on interdisciplinary work by engineers, mathematicians, and physicists. We derive a multiphysics model, that includes a simplification of the dominant physical processes, that is, transport of the reactive species in the gas phase and their deposition rates at the target layer.

The approximation between the mathematical parameters and physical parametes can be done by a regression method, such that we can verify the physical experiments. Such approximations help to study the physical experiments with simulation tools which are cheaper and can predict the more appropriate experiments which should be done to understand and control the physical processes.

In the following, we introduce the PE-CVD process and the important trends in modelling it.

A gas exposed to an electric field under low pressure conditions (<5Torr) results in a nonequilibrium plasma; see [3, 4]. Such ionized media, known as “cold” plasma or glow discharges, are powerful surface-modification tools in Material Science and Technology. Low-pressure plasmas allow to modify the surface chemistry and properties of materials compatible with a low-medium vacuum, through a PE-CVD process; see the applications in [4, 5].

Here, PE-CVD processes are attractive methods, because of their reproducible chemical processes that can be controlled by pressure, temperature, and by additional precursor gases. Such methods have been developed in recent years and there is an interest in producing high-temperature films; see [2].

We consider models that are related to mesoscopic scales [6], with respect to flows close to the wafer surface, where a wafer is the target material (e.g., metal or ceramic) for deposition, [2]. We assume that the wafer is a homogeneous medium and that the surface can be modelled as a porous medium [7].

The physical experiments are used to measure the influence of temperature, pressure and plasma power on deposition rates; see [8]. Here the plasma reactor chamber of an NIST GEC reference cell is used and for the hybrid ICP/CCP-RF plasma source, a double spiral antenna is used; see [8]. Such experiments are important but the variation of all possible parameters is very extensive.

Mathematically, we apply an interpolation between physical and mathematical parameters to verify our simulation model. Based on the smaller mathematical parameter space, we can allow much more experiments and obtain via the regression function the resulting parameters to the physical experiments. Such switching between numerical experiments and physical experiments reduces the experiments to a possible amount and we can optimize the deposition process.

The numerical results are discussed and applied to validation problems and real-life problems. We discuss an application involving the deposition of SiC thin films on a metallic plate.

The paper is organized as follows.

In Section 2, we present our mathematical model and a possible reduced model for further approximations. In Section 3, we discuss the physical experiments of the CVD process. The numerical methods of the transport-reaction equation and their parameter approximation to the physical model are described in Section 4. The numerical experiments are given in Section 5. In Section 6, we briefly summarize our results.

2. Mathematical Model

In the next, we discuss the derivation of our model.

We start by developing a multiphase model in the following steps:(i)standard transport model (one phase),(ii)flow model (flow field of plasma medium), (iii)multiphase model with mobile and immobile zones.

In each part of the model, we can refine the transport processes of the deposition gaseous species or reaction gaseous species depending on the influence of flow field, plasma zones and precursor gases.

A schematic test geometry of the CVD reactor is given in Figure 1.

Figure 1 

Vertically impinging CVD reactor.

2.1. Standard Transport Model

In the following, the models are discussed in terms of far-field and near-field problems, which take into account the scales of the model.

Two different types of model can be discussed:(1)convection-diffusion-reaction equations [6] (far-field problem), (2)Boltzmann-Lattice equations [9] (near-field problem).

The modelling is governed by the Knudsen Number, whereby the Knudsen number is a dimensionless number and defines the ratio of the molecular mean free path length to a representative physical length scale.𝜆Kn=𝐿, (2.1) where 𝜆 is the mean free path and 𝐿 is the representative physical length scale. This length scale could be, for example, the radius of a body in a fluid. Here, we deal with small Knudsen Numbers Kn≈0.01−1.0 for a convection-diffusion-reaction equation and a constant velocity field, whereas for, large Knudsen Numbers Kn≥1.0, we deal with a Boltzmann equation [2]. From modelling the gaseous transport of the deposition species, we consider the pure far-field model and assume a continuum flow field; see [10].

Such assumptions lead to transport equations that can be treated with a convection-diffusion-reaction equation owing to a constant velocity field, thus:𝜕𝜕𝑡𝑐+∇𝐹−𝑅_𝑔=𝑞(𝑥,𝑡),inΩ×[0,𝑡],𝐹=𝐯𝑐−𝐷∇𝑐,𝑐(𝑥,𝑡)=𝑐₀(𝑥),onΩ,𝑐(𝑥,𝑡)=𝑐₁(𝑥,𝑡),on𝜕Ω×[0,𝑡], (2.2)(2.3) where 𝑐 is the molar concentration of the reaction gases (called species) and 𝐹 is the flux of the species. 𝐯 is the flux velocity through the chamber and porous substrate [7]. 𝐷 is the diffusion matrix and 𝑅_𝑔 is the reaction term. The initial value is given as 𝑐₀ and we assume a Dirichlet boundary with the function 𝑐₁(𝑥,𝑡) sufficiently smooth. 𝑞(𝑥,𝑡) is a source function depending on time and space and represents the inflow of the species.

The parameters of the equation are derived as follows. Diffusion in the modified CVD process is given by Knudsen diffusion [11]. We consider the overall pressure in the reactor as 200 Pa and the substrate temperature (or wafer surface temperature) is about 600–900 K. The pore size in the homogeneous substrate is assumed to be 80 nm. The homogeneous substrate can be either a porous medium, for example, a ceramic material, see [11], or a dense plasma, assumed to be very dense and stationary; see [1]. For such media, we can derive diffusion based on Knudsen diffusion.

The diffusion is described by𝐷=2𝜖𝜇_𝐾𝜈𝑟,3𝑅𝑇 (2.4) where 𝜖 is the porosity, 𝜇_𝐾 is the shape factor of Knudsen diffusion, 𝑟 is the average pore radius, 𝑅 and 𝑇 are the gas constant and temperature, respectively, and 𝜈 is the mean molecular speed, given by𝜈=8𝑅𝑇𝜋𝑊, (2.5) where 𝑊 is the molar mass of the diffusive gas.

For homogeneous reactions during the CVD process, we consider a constant reaction of Si, Ti, and C species given as3Ti+Si+2C→Ti₃SiC₂, (2.6) where Si₃TiC₂ is a MAX-phase material, see [12], which deposits on the wafer surface. For simplicity, we do not consider the intermediate reaction with the precursor gases [1] and we assume we are dealing with a compound gas 3Ti+Si+2C; see [13]. Therefore we can concentrate on the transport of one species.

The reaction rate is then given by𝜆=𝑘_𝑟[3Si]^𝑀[Ti]^𝑁[2C]^𝑂[Si₃TiC₂]^𝐿, (2.7) where 𝑘_𝑟 is the apparent reaction constant and 𝐿,𝑀,𝑁,𝑂 are the reaction orders of the reactants.

A schematic overview of the one-phase model is presented in Figure 2. Here, the gas chamber of the CVD apparatus is shown, which is modelled by a homogeneous medium.

Figure 2 

Gas chamber of the CVD apparatus.

2.2. Flow Field

The flow field is derived in which the velocity is used for the transport of species. The velocity in the homogeneous substrate is modelled by a porous medium [14, 15]. We assume a stationary or low reactivity medium, for example, nonionized or low-ionized plasma or less reactive precursor gas. Further, the pressure can be assumed to have the Maxwell distribution [1]:𝑝=𝜌𝑏𝑇, (2.8) where 𝜌 is the density, 𝑏 is Boltzmann's constant, and 𝑇 is the temperature.

The model equations are based on mass and momentum conservation equations, where we assume conservation of energy. Because of the low temperature and low pressure environment, we assume the gaseous flow has a nearly liquid behaviour. Therefore, the derivation of the velocity can be given by Darcy's law:𝑘𝐯=−𝜇(∇𝑝−𝜌𝑔), (2.9) where 𝐯 is the velocity of the fluid, 𝑘 is the permeability tensor, 𝜇 is the dynamic viscosity, 𝑝 is the pressure, 𝑔 is the vector of gravity, and 𝜌 is the density of the fluid.

We use the continuum equation of the particle density and obtain the equation of the system, which is given as a flow equation:𝜕_𝑡(𝜙𝜌)+∇⋅(𝜌𝐯)=𝑄, (2.10) where 𝜌 is the unknown particle density, 𝜙 is the effective porosity, and 𝑄 is the source-term of the fluid. We assume a stationary fluid and consider only divergence-free velocity fields, that is,∇⋅𝐯(𝑥)=0,𝑥∈Ω. (2.11) The boundary conditions for the flow equation are given as𝑝=𝑝_𝑟(𝑡,𝛾),𝑡>0,𝛾∈𝜕Ω,𝐧⋅𝐯=𝑚_𝑓(𝑡,𝛾),𝑡>0,𝛾∈𝜕Ω, (2.12) where 𝐧 is the normal unit vector with respect to 𝜕Ω, where we assume that the pressure 𝑝_𝑟 and flow concentration 𝑚_𝑓 are prescribed by Dirichlet boundary conditions [15].

From the nearly stationary fluids, we assume that the conservation of momentum for velocity 𝐯 is given [15, 16]. Therefore, we can neglect the computation of the momentum for the velocity.

Remark 2.1. For flow through a gas chamber, for which we assume a homogeneous medium and nonreactive plasma, we have considered a constant flow [17]. A further simplification is given by the very small porous substrate, for which we can assume the underlying velocity in a first approximation as constant [2].

Remark 2.2. For a nonstationary media and reactive or ionized plasma, we have to take into account the relations for electrons in thermal equilibrium. Such a spatial variation can be considered by modelling electron drift. Such modelling of the ionized plasma is done with the Boltzman relation [1].

2.3. Multiphase Model: Mobile and Immobile Zones

More complicated processes such as retardation, adsorption, and dissipation processes, the gaseous species are modelled with multiphase equations. We take into account the fact that the concentration of species can be given in mobile and immobile versions, depending on their different reactive states; see [7]. From these behaviours, we have to model the transport and adsorbed state of species; see also Figure 3. Here, the mobile and immobile phases of the gas concentration are shown at a macroscopic scale in the porous medium.

Figure 3 

Mobile and immobile phases.

The model equations are given as combinations of transport and reaction equations (coupled partial and ordinary differential equations) as𝜙𝜕_𝑡𝑐^𝐿_𝑖+∇⋅(𝐯𝑐^𝐿_𝑖−𝐷^𝑒(𝑖)∇𝑐^𝐿_𝑖)=𝑔(−𝑐^𝐿_𝑖+𝑐^𝐿_{𝑖,𝑖𝑚})−𝜆_𝑖,𝑖𝜙𝑐^𝐿_𝑖+_{𝑘=𝑘(𝑖)}𝜆_𝑖,𝑘𝜙𝑐^𝐿_𝑘+𝑄_𝑖,𝜙𝜕_𝑡𝑐^𝐿_{𝑖,𝑖𝑚}=𝑔(𝑐^𝐿_𝑖−𝑐^𝐿_{𝑖,𝑖𝑚})−𝜆_𝑖,𝑖𝜙𝑐^𝐿_{𝑖,𝑖𝑚}+_{𝑘=𝑘(𝑖)}𝜆_𝑖,𝑘𝜙𝑐^𝐿_{𝑘,𝑖𝑚}+𝑄_{𝑖,𝑖𝑚}, (2.13)(2.14) where 𝑖=1,…,𝑀 and 𝑀 denotes the number of components.

The parameters in (2.13) are further described; see also [18].

The effective porosity is denoted by 𝜙 and describes the portion of the porosity of the aquifer that is filled with plasma, and we assume a nearly fluid phase. The transport term is indicated by the Darcy velocity 𝐯, which presents the flow direction and absolute value of the plasma flux. The velocity field is divergence-free. The decay constant of the 𝑖th species is denoted by 𝜆_𝑖. 𝑘(𝑖) therefore denotes the indices of the other species.

Remark 2.3. The concentrations in the mobile zones is modelled with convection-diffusion-reaction equations, see also Section 2.1, while the concentration in the immobile zones are modelled with reaction equations. These two phases represent the mobility of the gaseous species through homogeneous media, where the concentrations in the immobile zones are at least the lost amounts of depositable gases.

2.4. Simplified Model: Far-Field Model

We concentrate on a far-field model and assume continuum flow and that the transport equations can be treated by a convection-diffusion-reaction equation with a constant velocity field:𝜕𝜕𝑡𝑐+∇𝐹−𝑅_𝑔=0,inΩ×[0,𝑡],𝐹=𝐯𝑐−𝐷∇𝑐,𝑐(𝑥,𝑡)=𝑐₀(𝑥),onΩ,𝑐(𝑥,𝑡)=𝑐₁(𝑥,𝑡),on𝜕Ω×[0,𝑡], (2.15) where 𝑐 is the molar concentration and 𝐹 is the flux of the species. 𝐯 is the flux velocity through the chamber and porous substrate [7]. 𝐷 is the diffusion matrix and 𝑅_𝑔 is the reaction term. The initial value is given as 𝑐₀ and we assume a Dirichlet boundary with the function 𝑐₁(𝑥,𝑡) sufficiently smooth.

Remark 2.4. The focus on the dominant far field processes in the gas phase involving the reactive species reduces enormously the physical parameter space. Such a realistic reduction with respect to the experiments can also simplify the underlying mathematical model and concentrate on a defined number of experiments. Such experiments can validate the switching between physical and mathematical parameter space and allow to foresee the important processes in the gas phase.

3. Physical Experiments

The base of the experimental setup is the plasma reactor chamber of an NIST GEC reference cell. The spiral antenna of a hybrid ICP/CCP-RF plasma source was replaced by a double spiral antenna [8]. This reduces the asymmetry of the magnetic field due to the superposition of the induced fields of both antennas. Also, the power coupling to the plasma increases and enhances the efficiency of the source. A set of MKS mass flow-controllers allow any defined mixture of gaseous precursors. Even the flows of liquid precursors with high vapor pressure are controlled by this system. All other liquid and all solid precursors will be directly transported to the chamber by a controlled carrier gas flow. Besides the precursor flow, the density can also be changed by varying the pressure inside the recipient. Controlling the pressure is achieved with a valve between the recipient and vacuum pumps. In addition, a heated and insulated substrate holder was mounted. Thus, a temperature up to 800°C and a bias voltage can be applied to the substrate. While the pressure and RF power determine the undirected particle energy (plasma temperature), the bias voltage adds, only to the charged particles, energy directed at the substrate. Apart from the pressure and RF power control, the degree of ionization and number as well as size of molecular fractions can be controlled.

Altogether, this setup provides as free process parameters: (i)pressure (typically 10⁻¹-10⁻² mbar),(ii)precursor composition (TMS,TMS+H₂,TMS+O₂),(iii)precursor flow-rate (ranging from SCCM to SLM), (iv)RF Power (up to 1100 W),(v)substrate temperature (RT–800°C), (vi)bias voltage (DC, unipolar and bipolar pulsed, floating).

During all experiments, the process was observed by optical emission spectroscopy (OES) and mass spectroscopy (MS). The stoichiometry of the deposited films was analyzed ex situ in a scanning electron microscope (SEM) by energy dispersive X-ray analysis (EDX).

Realisation of Physical Experiments
The following parameters are used for the physical experiments. Such a reduction allows to concentrate on the important flow and transport processes in the gas phase. Further, we apply the underlying mathematical model (far field model, see Section 2.4) such that we can switch between physical and mathematical parameters.

Precursor: Tetramethylsilan (TMS),Substrate: VA-Steel,Film at substrate: SiC_𝑥(see Table 1).

Table 1 

We apply the parameters shown in Table 2 for interpolation of the substrate temperature.

Table 2 

For the substrate temperature and plasma power, see Table 3.

Table 3 

Remark 3.1. In the process, the temperature and power of the plasma are important and experiments show that these are significant parameters. Based on these parameters, we initialize the mathematical model and interpolate the flux and reaction parts.

4. Numerical Methods

In this section, we discuss the numerical methods. To accelerate our numerical methods, we combined the numerical and analytical parts in the solver processes.

4.1. Discretization and Solver Methods

For space-discretization of the PDEs, we apply finite-volume methods in mass conserved discretization schemes and for time-discretization of the resulting ODEs we apply Runge-Kutta methods or BDF methods. To accelerate the solver process, we combine the numerical and analytical parts of the solutions.

4.1.1. Discretization Method for Convection Equation

We deal with the following convection equation𝜕_𝑡𝑅𝑐−𝐯⋅∇𝑐=0, (4.1) where 𝑅 is the retardation factor and presents the retention of the concentration; see also (2.2). 𝐯 is the velocity. We have a simple boundary condition 𝑐=0 for the inflow and outflow boundary and the initial values are given as 𝑐(𝑥_𝑗,0)=𝑐⁰_𝑗(𝑥). We use a piecewise constant discretization method with the upwind discretization done in [19] and get𝑉_𝑗𝑅𝑐_𝑗^𝑛+1=𝑉_𝑗𝑅𝑐^𝑛_𝑗−𝜏^𝑛_{𝑘∈out(𝑗)}𝑣_𝑗𝑘𝑐^𝑛_𝑗+𝜏^𝑛_{𝑘∈in(𝑗)}𝑐^𝑛_𝑘𝑣_𝑘𝑗,𝑉_𝑗𝑅𝑐_𝑗^𝑛+1=𝑐^𝑛_𝑗(𝑅𝑉_𝑗−𝜏^𝑛𝜈_𝑗)+𝜏^𝑛_{𝑘∈in(𝑗)}𝑐^𝑛_𝑘𝑣_𝑘𝑗. (4.2) The explicit time discretization has to fulfill the discrete minimum-maximum property [19] and we get the following restriction on time steps:𝜏_𝑗=𝑅𝑉_𝑗𝜈_𝑗,𝜏^𝑛≤min_{𝑗=1,…,𝐼}𝜏_𝑗. (4.3) To obtain improved spatial discretization methods and apply larger time steps, we introduce a reconstruction with linear polynomials as a higher test-function in the next subsection.

4.1.2. Discretization Method for Convection-Reaction Equation Based on Embedded One-Dimensional Analytical Solutions

We apply Godunov's discretization method, compare for example, [20], and extend the formulation to the analytical solution of convection-reaction equations. We reduce the multidimensional equation to one-dimensional equations and solve each equation exactly. The one-dimensional solution is multiplied by the underlying volume to give the mass formulation. The one-dimensional mass is embedded in the multidimensional mass formulation and we obtain the discretization of the multidimensional equation.

The algorithm is given in the following manner:𝜕_𝑡𝑐_𝑙+∇⋅𝐯_𝑙𝑐_𝑙=−𝜆_𝑙𝑐_𝑙+𝜆_𝑙−1𝑐_𝑙−1,with𝑙=1,…,𝑚. (4.4) The velocity vector 𝐯 is divided by 𝑅_𝑙. The initial conditions are given by 𝑐⁰₁=𝑐₁(𝑥,0), else 𝑐⁰_𝑙=0for𝑙=2,…,𝑚 and the boundary conditions are trivial 𝑐_𝑙=0 for 𝑙=1,…,𝑚.

We first calculate the maximum time step for cell 𝑗 and concentration 𝑖 with the use of the total outward fluxes:𝜏_𝑖,𝑗=𝑉_𝑗𝑅_𝑖𝜈_𝑗,𝜈_𝑗=_{𝑘∈out(𝑗)}𝑣_𝑗𝑘. (4.5) We get the restricted time step with the local time steps of cells and their components𝜏^𝑛≤min_{𝑖=1,…,𝑚𝑗=1,…,𝐼}𝜏_𝑖,𝑗. (4.6) The velocity of the discrete equation is given by𝑣_𝑖,𝑗=1𝜏_𝑖,𝑗. (4.7) We calculate the analytical solution of the mass, compare for example, [18], and we get𝑚^𝑛_{𝑖,𝑗𝑘,out}=𝑚_𝑖,out(𝑎,𝑏,𝜏^𝑛,𝑣_1,𝑗,…,𝑣_𝑖,𝑗,𝑅₁,…,𝑅_𝑖,𝜆₁,…,𝜆_𝑖𝑚),_{𝑛𝑖,𝑗,rest}=𝑚^𝑛_𝑖,𝑗𝑓(𝜏^𝑛,𝑣_1,𝑗,…,𝑣_𝑖,𝑗,𝑅₁,…,𝑅_𝑖,𝜆₁,…,𝜆_𝑖), (4.8) where 𝑎=𝑉_𝑗𝑅_𝑖(𝑐^𝑛_{𝑖,𝑗𝑘}−𝑐^𝑛_{𝑖,𝑗𝑘′}), 𝑏=𝑉_𝑗𝑅_𝑖𝑐^𝑛_{𝑖,𝑗𝑘′}, and 𝑚^𝑛_𝑖,𝑗=𝑉_𝑗𝑅_𝑖𝑐^𝑛_𝑖,𝑗. Furthermore, 𝑐^𝑛_{𝑖,𝑗𝑘′} is the concentration at the inflow and 𝑐^𝑛_{𝑖,𝑗𝑘} is the concentration at the outflow boundary of cell 𝑗.

The discretization with the embedded analytical mass is calculated by:𝑚^𝑛+1_𝑖,𝑗−𝑚^𝑛_𝑖,rest=−_{𝑘∈out(𝑗)}𝑣_𝑗𝑘𝜈_𝑗𝑚_{𝑖,𝑗𝑘,out}+_{𝑙∈in(𝑗)}𝑣_𝑙𝑗𝜈_𝑙𝑚_{𝑖,𝑙𝑗,out}, (4.9) where 𝑣_𝑗𝑘/𝜈_𝑗 is the retransformation of the total mass 𝑚_{𝑖,𝑗𝑘,out} into the partial mass 𝑚_{𝑖,𝑗𝑘}. In the next time step, the mass is given as 𝑚^𝑛+1_𝑖,𝑗=𝑉_𝑗𝑐^𝑛+1_𝑖,𝑗 and in the old time step it is the rest mass for concentration 𝑖. The proof is given in [18]. In the next section, we derive an analytical solution for the benchmark problem, compare for example, [21, 22].

In the next subsection, we introduce the discretization of the diffusion-dispersion-equation.

4.1.3. Discretization of the Diffusion-Dispersion-Equation

We discretize the diffusion-dispersion-equation with implicit time-discretization and finite-volume method for the following equation:𝜕_𝑡𝑅𝑐−∇⋅(𝐷∇𝑐)=0, (4.10) where 𝑐=𝑐(𝑥,𝑡) with 𝑥∈Ω and 𝑡≥0. The diffusion-dispersion tensor 𝐷=𝐷(𝑥,𝐯) is given by the Scheidegger approach, compare for example, (4.15). The velocity is given as 𝐯. The retardation factor is 𝑅>0.0.

The boundary-values are denoted by 𝐧⋅𝐷∇𝑐(𝑥,𝑡)=0, where 𝑥∈Γ is the boundary Γ=𝜕Ω, compare for example, [23]. The initial conditions are given by 𝑐(𝑥,0)=𝑐₀(𝑥).

We integrate (4.10) over space and time and derive_Ω𝑗^𝑡𝑛+1_𝑡𝑛𝜕_𝑡𝑅(𝑐)𝑑𝑡𝑑𝑥=_Ω𝑗^𝑡𝑛+1_𝑡𝑛∇⋅(𝐷∇𝑐)𝑑𝑡𝑑𝑥. (4.11) The time-integration is done by the backward-Euler method and the diffusion-dispersion term is lumped, compare for example, [18]._Ω𝑗(𝑅(𝑐^𝑛+1)−𝑅(𝑐^𝑛))𝑑𝑥=𝜏^𝑛_Ω𝑗∇⋅(𝐷∇𝑐^𝑛+1)𝑑𝑥. (4.12) Equation (4.12) is discretized over space by using Green's formula._Ω𝑗(𝑅(𝑐^𝑛+1)−𝑅(𝑐^𝑛))𝑑𝑥=𝜏^𝑛_Γ𝑗𝐷𝐧⋅∇𝑐^𝑛+1𝑑𝛾, (4.13) where Γ_𝑗 is the boundary of the finite-volume cell Ω_𝑗. We use an approximation in space, compare for example, [18].

The spatial integration of (4.13) is done by the mid-point rule over finite boundaries and is given by𝑉_𝑗𝑅(𝑐_𝑗^𝑛+1)−𝑉_𝑗𝑅(𝑐^𝑛_𝑗)=𝜏^𝑛_{𝑒∈Λ𝑗𝑘∈Λ𝑒𝑗}|Γ^𝑒_𝑗𝑘|𝐧^𝑒_𝑗𝑘⋅𝐷^𝑒_𝑗𝑘∇𝑐^𝑒,𝑛+1_𝑗𝑘, (4.14) where |Γ^𝑒_𝑗𝑘| is the length of the boundary-element Γ^𝑒_𝑗𝑘. The gradients are calculated with the piecewise finite-element function 𝜙_𝑙, see [24], and we obtain∇𝑐^𝑒,𝑛+1_𝑗𝑘=_{𝑙∈Λ𝑒}𝑐_𝑙^𝑛+1∇𝜙_𝑙(𝐱^𝑒_𝑗𝑘). (4.15) We use the difference notation for the neighbour points 𝑗 and 𝑙, compare for example, [25], and obtain the discretized equation:𝑉_𝑗𝑅(𝑐_𝑗^𝑛+1)−𝑉_𝑗𝑅(𝑐^𝑛_𝑗)=𝜏^𝑛_{𝑒∈Λ𝑗𝑙∈Λ𝑒⧵{𝑗}}(_{𝑘∈Λ𝑒𝑗}|Γ^𝑒_𝑗𝑘|𝐧^𝑒_𝑗𝑘⋅𝐷^𝑒_𝑗𝑘∇𝜙_𝑙(𝐱^𝑒_𝑗𝑘))(𝑐_𝑗^𝑛+1−𝑐_𝑙^𝑛+1), (4.16) where 𝑗=1,…,𝑚.

4.2. Interpolation and Regression of Experimental Dates

To simulate the physical experiments with the assumed model, we have to approximate the parameters of the numerical model. We apply interpolation and regression schemes to approximate between the mathematical and physical parameters.

Here, we concentrate on the reaction rates of species Si, C, and H.

The physical data of temperature and pressure are used and validation simulations are done to obtain the rate of deposition.

Next we have to interpolate the parameters of the numerical model.

(1) Lagrangian Interpolation
We assume an interpolation at Ω=[𝑎₁,𝑏₁]×⋯×[𝑎_𝑑,𝑏_𝑑]. 𝑇=_𝜈∈𝐾𝑓(𝑥^𝑡_𝜈)𝐿^𝑡_𝜈, (4.17) where the Lagrangian function is given as 𝐿^𝑡_𝜈(𝑥)=𝜋^𝑑_𝑖=1𝜋^𝑚_{𝜇=0,𝜇≠𝜈𝑖}𝑥_𝑖−𝑥^{[𝑎𝑖,𝑏𝑖]}_𝜇𝑥_𝜈𝑖𝑖^]_𝑖^[𝑎,𝑏−𝑥_{𝜇𝑖𝑖][𝑎,𝑏}. (4.18)

(2) Linear Regression (Least squares approximation)
Here, we have points with values and we assume we have the best approximation by minimization: 𝑆=^𝑚_𝑘=1(𝑦_𝑘−𝐿_𝑛(𝑥_𝑘))², (4.19) where 𝑚≥𝑛 and 𝐿_𝑛 is a function that is constructed by the least squares algorithm; see [26].

Remark 4.1. To apply larger parameter spaces, we can generalise to multivariate regression methods see [27]. Here we compute approximations between higher dimensional matrix spaces.

5. Numerical Experiments

For all the experiments, we have the following parameters of the model, discretization and solver methods.

We apply interpolation and regression methods to couple the physical parameters to the mathematical parameters; see Figure 4 and Table 4.

Table 4 

Physical and mathematical parameters.

Figure 4 

Coupling of physical and mathematical parameter space.

Parameters of the Equation
In the following, we list the parameters for our simulation tool UG; see [28]. The software toolbox has a flexible user interface to allow a large number of numerical experiments and approximations to the known physical parameters. The model parameters are given in Table 5.

Table 5 

Model-Parameters.

Discretization Method
Finite-volume method of 2nd order are given with theparameters in Table 6.

Table 6 

Time discretization parameters.

Time Discretization Method
Crank-Nicolson method of 2nd order are given with theparameters in Table 7.

Table 7 

Spatial discretization parameters.

Solver Method
In the following, we deal with test examples which are approximations of physical experiments. The parameters of the solvers are given in Table 8.

Table 8 

Solver methods and their parameters.

5.1. Test Experiment 1: Interpolation with Temperature

In the test example, we deal with the following reaction:2SiC+4H→_𝜆SiC+CH4+Si. (5.1) Here, we have physical experiments and approximate temperature parameters 𝑇=400,600,800.

We compute the SiC∶C ratio at a given temperature 𝑇=400,600,800 with the UG program and fit to the parameter 𝜆.

We use a Lagrangian formula to compute 𝜆 for the new temperatures 𝑇=500,700 and apply the ratio of the simulated new parameters. These values can be given back to physical experiments; see Table 9.

Table 9 

Computed and experimental fitted parameters from UG simulations.

One Source
The parameters of the source concentration are givenin Table 10. In Figure 5, we present the concentration of one point source at (50,20) with number of time steps being equal to 25.

Table 10 

Parameters of source concentration.

(a)

(b)

Figure 5 

One point source at (50,20) with number of time steps being equal to 25.

In Figure 6, we show the deposition rates of one point source at (50,20), with number of time steps being equal to 25. The rate of concentrations is given in Table 12.

Figure 6 

Deposition rates in case of one point source at (50,20) with number of time steps being equal to 25.

Nine Point Sources
In this experiment, we apply nine point sources.
The parameters of the source concentration are givenin Table 11. In Figure 7, we present the concentration with nine point sources in a short time.
In Figure 8, we show the deposition rate with the nine point sources, with number of time steps being equal to 25. The rate of concentrations is given in Table 13.

Table 11 

Parameters of source concentration.

Table 12 

Rate of concentration.

Table 13 

Rate of concentration.

(a)

(b)

Figure 7 

Nine point sources with number of time steps being equal to 25.

Figure 8 

Deposition rates in case of nine point sources, with number of time steps being equal to 25.

81-Point Sources
In this experiment, we apply 81 point sources.
The parameters of the source concentration are givenin Table 14. In Figure 9, we present the concentration with 81 point sources.
In Figure 10, we show the deposition rates with 81 point sources, with number of time steps equal to 80. The rate of concentrations is given in Table 15.

Table 14 

Parameters of source concentration.

Table 15 

Rate of concentration.

(a)

(b)

Figure 9 

81 point sources with number of time steps equal to 80.

Figure 10 

Deposition rates in case of 81 point sources, with number of time steps equal to 80.

Line Source
In this part, we will perform experiments with a line source. The parameters of the line source are given in Table 16.
In Figure 11, we present the results with a line source, 𝑥∈[5,95],𝑦∈[20,25] with number of time steps equal to 25.
In Figure 12, we show the deposition rates with (a) the line source, 𝑥 is from 5 to 95, and 𝑦 is from 20 to 25. The rate of concentrations is given in Table 17.

Table 16 

Parameters the source concentration.

Table 17 

Rate of concentration.

(a)

(b)

Figure 11 

Line source with number of time steps being equal to 25.

Figure 12 

Deposition rates in case of line source, $x\in [\mathrm{5,95}],y\in [\mathrm{20,25}]$.

5.2. Test Experiment 2: Interpolation with Temperature and Power

In the next experiment, we fit the mathematical parameters to the temperature and power of the physical experiments.

We deal with the reaction:2SiC+4H→_𝜆SiC+CH4+Si. (5.2) In this case, we have a table which has values for the temperature and power of the plasma and for the ratio between the sources.

We have to interpolate 𝜆 to the physical parameters temperature 𝑇 and power of plasma 𝑃. In Table 18, the interpolated parameters are given.

Table 18 

Computed (C) and experimentally fitted (F) parameters from UG simulations.

One Source
The parameters of the source concentration are givenin Table 19. In Figure 13, we present the concentration of one point source at (50,20) with number of time steps being equal to 25. In Figure 14, we show the deposition rates with one point source at (50,20), with number of time steps being equal to 25. The rate of concentrations is given in Table 20.

Table 19 

Parameters of source conce𝑛𝑡𝑟𝑎𝑡𝑖𝑜𝑛.

Table 20 

Rate of concentration.

(a)

(b)

(c)

(d)

Figure 13 

One point source at (50,20) with number of timesteps being equal to 25.

Figure 14 

Deposition rates in case of one point source at (50,20), with number of time steps being equal to 25.

Nine Point Sources
In this experiment, we apply nine point sources.
The parameters of the source concentration are given in Table 21. In Figure 15, we present the concentration of the nine point sources in a short time.
In Figure 16, we show the deposition rates of nine point sources, with number of time steps being equal to 25.

Table 21 

Parameters of source concentration.

(a)

(b)

(c)

(d)

Figure 15 

Nine point sources with number of time steps being equal to 25.

Figure 16 

Deposition rates in case of nine point sources, with number of time steps being equal to 25.

81-Point Sources
In this experiment, we apply 81 point sources.
The parameters of the source concentration are givenin Table 23. In Figure 17, we present the concentration with 81 point sources.
In Figure 18, we show the deposition rates with 81 point sources, with number of time steps equal to 100. The rate of concentrations is given in Table 22.

Table 22 

Rate of concentration.

Table 23 

Parameters of source concentration.

(a)

(b)

(c)

(d)

Figure 17 

81 point sources with number of time steps equal to 100.

Figure 18 

Deposition rates in case of 81 point sources, with number of time steps equal to 100.

Line Source
In this part, we will perform an experiments with a line source, 𝑥∈[5,95],𝑦∈[20,25].
In Figure 19, we present the result of the line source, 𝑥∈[5,95],𝑦∈[20,25], with number of time being steps equal to 30.
In Figure 20, we show the deposition rates with a line source, with 𝑥 from 5 to 95, and 𝑦 from 20 to 30.

(a)

(b)

(c)

(d)

Figure 19 

Line source with number of time steps being equal to 25.

Figure 20 

Deposition rates in case of line source, $x\in [\mathrm{5,95}],y\in [\mathrm{20,25}]$.

5.3. Test Experiment 3: Regression with Temperature and Power

In the next experiment, we apply a more flexible approximation method to obtain the parameters of the mathematical method. We apply a regression and fit to all the physical parameters, because we are not restricted to a given interpolation grid.

The reaction is given as

𝐴→𝐵 and 𝐵→C and we apply to 2SiC+4H→SiC+CH4+Si.

We computed the ratio SiC∶C for temperatures 𝑇=400,600,800 and plasma powers 100,300,500,900 and fit the simulated ratio given by the UG program to the mathematical model with a reaction parameter 𝜆.

We use linear regression, see Section 4, and compute 𝜆 for the new temperatures 𝑇=450,500,800 and apply the ratio of the simulated new parameters. These values can be given back to the physical experiments; see Table 24.

Table 24 

Rate of concentration.

Table 25 

Parameters of source concentration.

Table 26 

Rate of concentration.

Table 27 

Parameters of source concentration.

Table 28 

Parameters of source concentration.

Table 29 

Rate of concentration.

Table 30 

Parameters of source concentration.

Table 31 

Rate of concentration.

Table 32 

Parameters of source concentration.

Table 33 

Rate of concentration.

One Source
Here, we take points sources.
In Figure 21, we present the concentration in the one-point source experiment.
In Figure 22, we show the deposition rates in the one-point source experiment.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 21 

One point source experiment.

Figure 22 

Deposition rates in case of one-point source experiment.

Nine Point Sources
Here, we take nine point sources of both concentrations.
In Figure 23, we present the concentration of the nine point sources experiment.
In Figure 24, we show the deposition rates of the nine point sources experiment.

(a)

(b)

(c)

Figure 23 

Nine-point source experiment.

Figure 24 

Deposition rates in case of nine-point source experiment.

81-Point Sources, 𝑡_half=2⋅10⁸
In this first experiment, the value of temperature is 400 deg C and 𝜆 is 0.5⋅10⁻⁸. Here, we take the concentration of SiC as a point source, and the concentration of H as a line source.
In Figure 25, we present the concentration of the 81-point source experiment.
In Figure 26, we show the deposition rates in the 81 point source experiment.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 25 

81-point source experiment.

Figure 26 

Deposition rates in case of 81 point source experiment.

Remark 5.1. The regression method is more flexible for approximating to the physical parameters. We obtain numerical results for different parameter studies, that are fitted to the physical experiments. The first test examples with multiple sources and temperature regions which are interesting to physicists are simulated. Here, we have coupled a mathematical model to a physical experiment and studied the near region of the deposition process.

6. Conclusions

We present a numerical simulation of a CVD process for depositing SiC films. Based on the different scales of physical and mathematical experiments, we apply a parameter approximation to fit the physical experiment to the mathematical experiment. Numerical approximations to the experimental data include the new parameters of the mathematical model. Such experiments allow to reduce the number of physical experiments to an acceptable number and give engineers and experimentalists a mathematical tool for predicting complex physical processes.

The first numerical results show predictions of physical experiments with a transport-reaction equation of the deposition process.

The temperature of the target and power of the plasma are chosen in such manner that the simulation results can help find the optimal deposition conditions. Furthermore, multiple sources obtain the best results with homogeneous layer deposition.

Such numerical simulations help to predict the deposition rates of the underlying film, for example, SiC. In future, we will analyze the validity of the models with more complicated precursor gases. Here the outstanding multivariate analysis will be important for approximating a large number of parameters.

Nomenclature

𝜙:	effective porosity [−]
𝑐𝐿𝑖:	concentration of 𝑖th gaseous species in plasma chamber
𝑐𝐿𝑖,𝑖𝑚:	concentration of 𝑖th gaseous species in immobile zones of plasma chamber phase [mol/mm3]
𝐯:	velocity in plasma chamber [mm/nsec]
𝐷𝑒(𝑖):	element-specific diffusion-dispersion tensor [mm2/nsec]
𝜆𝑖,𝑖:	decay constant of 𝑖th species [1/nsec]
𝑄𝑖:	source term of 𝑖th species [mol/(mm3 nsec)]
𝑔:	exchange rate between mobile and immobile concentrations [1/nsec].