PHOENICS-CVD is a special-purpose code developed to simulate chemical vapour deposition (CVD) equipment, particularly in the semiconductor industry. The technical features of the code are described; these relate primarily to transport mechanisms, gas-phase and surface chemistry, surface-to-surface radiation and plasma effects. Implementation is by means of input and data files, together with a user-friendly menu system in which much of the set-up is automated . Extensive testing and validation of the code has been carried out.

Chemical Vapour Deposition (CVD) is a technique that is widely used, particularly in the semiconductor industry. There is a continuous demand for scale-up of reactors to produce larger wafers; at the same time the feature sizes are reducing, bringing ever increasing demands on film quality and process performance. As a result, the life cycle of CVD equipment is short and new equipment has to be developed every few years. Conventionally, such requirements have been met by empirical trial-and-error methods, straight-forward scale-up of existing equipment and ad hoc modifications to existing processes. However, computer simulation offers the possibility of a more fundamental understanding of the essential physics and chemistry, leading to a systematic approach to reactor and process design.

PHOENICS-CVD was developed to meet this challenge. It was created by a consortium of organisations within the ACCESS-CVD Project, with financial support from the European Commission through its ESPRIT programme. The partners in the Project were CHAM (co-ordinator), Technical University of Delft, Fraunhofer Institut fuer Integrierte Schaltungen (IIS-B), Siemens Corporate Research and Development and ASM International.

The objectives of the ACCESS-CVD Project were:

This document concentrates on the first of these objectives; the models that were developed for inclusion in the code are described in the following Section. For additional information the reader is referred to the special December 1995 issue of The PHOENICS Journal, which was exclusively devoted to the ACCESS-CVD Project; further references to original source material can be found there.

Various basic assumptions are made about the nature of the gas flow in the reactor. The first is that the gas mixture will behave as a continuum, i.e. that the mean free path for the molecules is much smaller than the characteristic dimensions of the reactor; this effectively imposes a lower limit on the operating pressure of, typically, about 30Pa. The assumption is made that ideal gas behaviour is a suitable approximation; the code has also been designed wi th laminar flow in mind, although this is not a requirement. Viscous heating due to dissipation is also neglected. All these limitations are entirely compatible with most CVD applications, although some very low pressure reactors should not be simulated using PHOENICS-CVD.

The fundamental equations to be solved are those representing conservation of mass, momentum and energy, coupled with a similar equation for each of the individual gas species. The mass and momentum equations take the conventional CFD form; the energy and species equations have additional contributions to take account of chemical reactions and diffusion of species within the gas mixture.

Species concentrations are expressed as dimensionless mass fractions, w _i; diffusion is expressed in terms of diffusive mass fluxes

(1)

where v_i is the i^th species velocity and v is the mass averaged mixture velocity. Diffusion arises from concentration gradients (ordinary diffusion) and also from temperature gradients (thermal or Soret diffusion); j_i therefore contains two separate contributions.

The conservation equation for the i^th species is then

(2)

where the terms on the right hand side represent, respectively, convection, diffusion and gas phase chemistry; m_i is molecular mass, is a stoichiometric coefficient and summation is over the reactions, with forward and backward reaction rates given, in molecular terms, by and .

The energy equation contains the standard convection and diffusion terms, augmented by the extra contributions referred to above:

(3)

where H_i is the species molar enthalpy and c_p is the mixture specific heat. The third term on the right hand side represents diffusion of energy associated with the inter-diffusion of species and the final term the energy associated with the creation and destruction of species by the chemical reactions.

Ordinary diffusion fluxes in a binary mixture are given by Fick’s Law:

(4)

where D₁₂ is the binary diffusion coefficient of the gas pair. The same approach can be applied to a multicomponent mixture in which a single carrier species dominates.

For non-dilute mixtures this simple formulation can no longer be utilised; instead, the diffusive fluxes of each species is related to all the concentration gradients through the Stefan-Maxwell equations:

(5)

where M is the average molar mass of the mixture. For implementation in PHOENICS-CVD it is more convenient to rearrange these equations to obtain expressions for each diffusion flux in terms of the others; this retains a conventional diffusion term in the species equations (useful for stability) with other terms being added as sources. A further advantage of the Stefan-Maxwell approach is that binary diffusion coefficients are used: these are simple functions of pressure and temperature; alternative methods require multicomponent diffusion coefficients, which are more complicated functions and correspondingly more expensive to evaluate.

As an alternative to the Stefan-Maxwell equations an approximate expression for multicomponent diffusion fluxes has been developed by Wilke; it retains the form of Fick’s Law,

(6)

but with an effective diffusion coefficient:

(7)

where x represents molar fraction. The Wilke approximation is accurate for dilute mixtures and offers computational savings when compared with Stefan-Maxwell.

The philosophy within PHOENICS-CVD has been to offer a range of alternatives to the user where these can save significant computational time. Fick, Wilke and Stefan-Maxwell formulations are therefore all included in the code.

Thermal diffusion causes small, light molecules to diffuse towards hot regions while larger, heavier molecules move the other way. This effect can be important in cold-wall reactors which may have large temperature gradients. The mass fluxes are given by

(8)

where the thermal diffusion coefficient is a function of gas composition and temperature, but independent of pressure.

For each gas species such properties as specific heat, heat of formation and standard entropy are needed in the energy equation given above (3); mixture specific heat is also required. These properties are functions of temperature and are usually defined by polynomial coefficients. PHOENICS-CVD follows the conventions adopted in the Chemkin thermodynamics database: for each species two sets of seven coefficients are required, one for a low temperature range and one for a high range.

Mixture specific heat is simply defined by mass averaging.

Solution of the conservation equations requires mixture properties such as density, viscosity and thermal conductivity to be known as functions of local conditions and composition.

Experimental data are available for only a few common gases and these properties are derived from kinetic theory for the individual gas species; Lennard-Jones parameters are required and these are taken from a data file whose format is again compatible with the corresponding Chemkin database.

Density poses no problems and is derived straightforwardly from the ideal gas law. Viscosity of a single gas is given by

(9)

where N_A is Avogadro’s number and is a tabulated function depending on the Lennard-Jones parameters for the species in question; the mixture viscosity is calculated from

(10)

with

. (11)

Thermal conductivity is given by

(12)

for monatomic gases; the Eucken correction is used for polyatomic gases:

. (13)

Mixture conductivity is derived from expressions having the same form as those used for viscosity.

Binary diffusion coefficients for the gas pairs are also obtained from kinetic theory, although a correction factor suggested by Wilke and Lee is also included; the result is

(14)

where is the average collision diameter for the pair and _Dij(T) is a tabulated function of the Lennard-Jones parameters and temperature.

Thermal diffusion modelling requires the use of multicomponent diffusion coefficients that are dependent on mixture composition. Derivation from kinetic theory is computationally expensive and increases as the cube of the number of gas species present; the problem lies in the evaluation of large determinants having elements which are themselves complicated functions of gas composition and temperature. As an alternative to this, a simpler model proposed by Clark Jones has also been inc luded in PHOENICS-CVD; use of this method results in

(15)

where the thermal diffusion factors (a_ij) are functions of molecular masses and Lennard-Jones parameters.

The user can select either the exact or the Clark Jones approach; as a further computational economy each formulation can be simplified by the use of the rigid elastic spheres approximation instead of the Lennard-Jones potential.

PHOENICS-CVD is equipped with two data files for transport and thermodynamic parameters. The user is able to add new species with complete freedom, following a standard format.

The fundamental process in Chemical Vapour Deposition is the chemistry that brings about the creation of the required film on the solid surface; inevitably this involves surface chemistry but generally there are also gas-phase reactions. Many simulation studies have relied on very simple models for the chemistry, with rate constants estimated from experimental data; while this approach might yield accurate predictions for the process conditions under which th e data were collected it is of limited value in the main field of predictive modelling. A deeper understanding of chemistry is also essential for studying conformality, dopant incorporation and the important effects that intermediates can have on process characteristics such as selectivity and deposit morphology. More detailed chemistry models have appeared recently, but they have typically only been applied in conjunction with idealised reactor configurations.

PHOENICS-CVD has been designed to facilitate the handling and solution of more sophisticated chemistry formulations, in addition to the simpler models, for both gas-phase and surface reactions. Numerical modelling can then be used to gain fundamental insight into the underlying phenomena.

Although not contributing directly to the deposition process itself, gas-phase chemistry plays an important role, particularly in atmospheric pressure reactors; the consumption of reactants and the formation of intermediaries may easily influence film quality.

The k^th reversible gas-phase reaction is written as

(16)

where ||n _ik|| = 0.5(|n _ik| +n _ik ).

In modelling terms the gas-phase chemistry contributes the final terms to the species and energy equations ((2) and (3)). The net reaction rate, - is obtained from

(17)

where k_k and k_-k are the forward and reverse reaction rate constants respectively. These constants are strongly dependent on temperature and independent of pressure at sufficiently high pressures. At lower pressures though there is pressure dependence.

A number of alternative formulations for the forward reaction rate constant have been included in PHOENICS-CVD. The simplest is the extended Arrhenius expression:

k_k = A Tb exp(-E_A/RT) p^c . (18)

More complicated alternatives are the Lindemann form which blends two Arrhenius expressions, one for low and one for high pressure. Further complexity is also available in the Troe form, based on nine or ten parameters.

For reversible reactions the reverse reaction rate constant is calculated from equilibrium thermochemistry, utilising the equilibrium constant for the reaction; this is derived from polynomial fits for standard heat of formation, standard entropy and specific heat.

Surface chemistry can be written in the same way as the gas-phase reactions described above:

(19)

where the A_i (i=1,N) are the gaseous reactants and products and the B_j (j=1,M) are the bulk, or adsorbed reactants and products; s _il and c _jl are the stoichiometric coefficients for the l^th surface reaction.

The growth rate is simply calculated from the reactions that deposit the bulk species:

(20)

where m_s is the molar mass of the bulk species and r _s is its density

Surface reactions did not appear explicitly in the species and energy equations ((2) and (3)) earlier. This is because they are included as boundary conditions rather than source terms.

2.2.2.1 Semi-empirical surface chemistry models

In many chemistry models the surface chemistry is simplified into one or more irreversible reactions with only gaseous reactants and, apart from the deposited species, only gaseous products. Equation (19) then simplifies to

. (21)

Semi-empirical methods have been developed for the definition of the surface reaction rates and two of these have been included as built-in options. The first involves the ‘reactive sticking coefficient’, defined as the fraction of gas molecules colliding with the surface that contribute to the deposited film. The reaction rate then becomes

(22)

where p_A is the partial pressure of the reactant and T_s the temperature at the surface; the sticking coefficient is g _A (£ 1) and the remainder of the expression comes from kinetic theory. Arrhenius expressions can be used for the sticking coefficient, specified via the chemistry data file.

A second built-in formulation is that by Langmuir-Hinshelwood; here the reaction rate is in the general form

. (23)

Again, the coefficients are specified in the data file.

2.2.2.2 Detailed surface chemistry models

The simplified surface chemistry models described above take no account of the adsorption of reactants and intermediaries at the surface, their decomposition and the desorption of reaction products. Although several more sophisticated models have been suggested in recent years there is still no standard form for them. In general, therefore, the inclusion of new models requires user-coding to be added to PHOENICS-CVD; this facility is provided and a number of specific models have been incl uded in the code to help users.

The solution strategy for detailed surface chemistry models is rather different from that usually employed. At the surface the fractional coverages of the adsorbed species are calculated using a separate solver for the simultaneous calculation of stiff surface chemistry; the resultant fluxes of reactants, intermediaries and products are then imposed in the normal way as boundary conditions to the transport equations for the gas species.

The solution of the transport equations for CVD processes is often hampered by numerical difficulties. These result from the chemistry source terms which lead to strong coupling between the mass fractions of different species at a single location; transport equations are usually solved on the basis that the strongest links are given by the convection and diffusion terms which link concentrations throughout the solution domain. Additionally, a single chemistry model is likely to contai n reactions with vastly differing time scales; this ‘stiffness’ can impose unacceptable limitations on the convergence that can be achieved.

Various techniques have been suggested to alleviate these problems, but they cannot usually be applied to two- and three-dimensional elliptic problems. In PHOENICS-CVD convergence is improved by the linearisation of the chemistry source terms which increases the stability of the solution method. The stiffness is addressed by the provision of an automatic under-relaxation based on local creation and destruction rates. In extreme cases this may still not be sufficient and an alternative solver has been provided; this solves the species equations on a point-by-point basis, using a Newton-Raphson solver for the simultaneous solution of all the species in a single computational cell. Such a technique is well-suited for handling large sets of chemical reactions with widely varying time scales; however, because it is a ‘local’ method convergence may still be slow because the convective and diffusive contributions are less well dealt with. In some cases the use of both solvers in turn may be beneficial.

Specification of the chemistry to be adopted is made easy by the use of a data file. For each of the numbered reactions the user is able to specify the chemical reaction (in terms of species and stoichiometric coefficients), the nature of the model to be adopted (e.g. kinetic limited, diffusion limited, Langmuir-Hinshelwood) and the appropriate values of the parameters in the model.

Radiation is a major element in heat transfer within CVD reactors, in addition to convection, conduction and diffusion. The nature of radiation is such that it does not fit easily into the framework of computational fluid dynamics, because of the links it introduces between remote parts of the solution domain. Radiation terms therefore appear as source terms calculated in a totally different manner from the other mechanisms of heat transfer. CVD gases are typ ically transparent to radiation and a surface-to surface model is therefore sufficient. There are two basic formulations: viewfactor based or Monte Carlo. Both approaches were successfully adopted within the ACCESS-CVD Project, although the Monte Carlo coding was restricted to axisymmetric geometries. Time constraints made it impossible for the Monte Carlo coding to be included in the first commercial release of PHOENICS-CVD; it is intended to add it to a later code update.

In each approach the surfaces in the reactor are divided into zones within which the temperature and surface optical properties are assumed to be constant. The radiative heat source to the i^th zone is given by

(24)

where T_j is the surface temperature of the j^th zone, s is the Stefan-Boltzmann constant and R is the radiative exchange matrix. This heat source is used together with the other fluxes to calculate the zone surface temperature iteratively as the general solution proceeds. It is the radiative exchange matrix that is calculated by the viewfactor or Monte Carlo approach.

The viewfactor from the j^th to the i^th thermal zone is given by

(25)

where r is the distance between the points in each zone, n_i and n_j are the two surface normals and A_j is the surface area of the j^th zone; is the transmittance from j to i in the l^th spectral band. If all surfaces are opaque the viewfactors are purely geometrical properties of the reactor; the introduction of semitransparent solids means that the viewfa ctors must be recalculated if the degree of transparency changes, e.g. as a result of a changing surface temperature.

Calculation of the energy flux from j to i in the l^th spectral band is calculated using the Gebhard factors

(26)

which take account of direct and multiply reflected radiation; r^(l)_m and a^(l)_m are the reflectivity and absorptivity of the m^th zone in the l^th band. The net flux in the l^th band is then given by

(27)

where is the fraction of black body radiation in the l^th band. The total radiative heat source, as given above, includes the full radiation exchange matrix which is the sum of the matrices for each spectral band.

The surface temperature of each zone is calculated iteratively as the solution proceeds and the surface temperature used for the viewfactors should take account of this. For simplicity a user-specified value for each zone is used to avoid frequent recalculation of the radiation exchange matrix. Greater accuracy can be achieved if the specified surface temperatures are modified manually and the solution restarted; typically a single iteration of this type is sufficient.

The calculation of the viewfactors between zones is carried out by means of a ray-tracing algorithm which detects obstructions (opaque or semitransparent). The geometrical viewfactors are then derived using a modified double integral method which takes account of singular integrals when zones are in contact. For computational efficiency the spectral variation of surface properties is implemented in a banded manner; this is sufficient for most purposes.

The Monte Carlo method is the same in principle; as before a radiative exchange matrix is generated from which heat sources can be determined. However, in this case the calculation is based directly on rays traced from each zone in turn until, after reacting with other surfaces, they reach their final destinations. A large number of rays are used, sampled from the appropriate angular and spectral distributions for the emitting (and any intermediate reflecting and transmitting) surface s, so that a statistically significant picture can be built up, showing how much energy from a surface will be received at each other zone. An advantage of the Monte Carlo approach is that it permits specular reflection to be simulated; this can be particularly useful for the highly reflective surfaces in reactor lamphouses. Furthermore, the spectral variation of properties can be exact, rather than in discrete bands.

A data file is used to specify the optical properties of any materials that are used in the simulation. The file contains, for each numbered material, a set of coefficients used in the calculation of the real and imaginary parts of the refractive index for each of sixty spectral bands; these cover the wavelength range from 10^-7m to 10^-4m, representing temperatures from 250K to 6000K. Simplified non-spectral modelling is also permitted and for this the data file c ontains constant values of emissivity and reflectivity. Additional materials can easily be added by users.

In recent years the use of plasma-enhanced chemical vapour deposition (PECVD) has become more common; the ability to provide the energy required for CVD chemistry without the undesirable side-effects of high thermal stress and poor controllability that are often a problem in conventional, thermally-energised systems makes it an attractive option for a range of processes.

The modelling of plasma-based chemistry is inherently more complex than the modelling of comparable neutral gas processes. Apart from the need to model additional variables associated with the electrons, ions and electric field, there is the notably non-equilibrium state of the plasma: the energy is taken up much more readily by the light electrons than the heavier ions and neutrals which typically remain close to the ambient temperature.

Various methods for plasma modelling exist, usually based either on fundamental laws of physics or engineering approaches containing mainly empirical information. The first category are typically very computer intensive while the second yield only qualitative results. The perceived need within the ACCESS-CVD Project was for a model that was simple enough to be compatible with the computational times for the rest of the simulation but still able to produce valuable quantitative results. Th e outcome was the Eddy Drift-Diffusion Model (EDDM) which has been implemented in PHOENICS-CVD. This model is explicitly geared towards capacitively coupled discharges in the PECVD regime; within those limits it is equivalent to the more fundamental fluid-dynamical approach, from which it is derived, but without such demanding computational requirements.

Unlike general kinetic models which describe a plasma in terms of distribution functions, fluid-dynamical plasma models use macroscopic variables such as densities, bulk velocities and pressures. Such a model for an RF discharge would require equations for electron density and energy (or temperature), density and velocity for each type of ion and also electrical potential. With appropriate boundary conditions these partial differential equations would describe the evolution of a low p ressure discharge under an applied electric field. For periodic RF excitation it would be necessary to solve the transient equations until the asymptotic limit cycle was reached; the transient effects, typically lasting only a few milliseconds, are of little interest for CVD plasma processes.

This approach has been shown to be very reasonable, but extremely slow. One reason for this is that stability and accuracy require time steps in the transient calculation which are small compared with the driving RF frequency; as a large number of cycles are required before the limit cycle is reached this leads to unreasonable run times. This is another example of numerical stiffness (see chemistry section above): the electrons and the ions are described within the same framework, yet the dynamical phenomena act on vastly different time scales.

The Effective Drift-Diffusion Model overcomes the stiffness problem by making use of additional simplifications to the fluid-dynamical equations; further, the limit cycle is calculated directly, with no attempt to model the transient evolution. The simplifications are based on the very circumstances which have caused the numerical stiffness in the first place: the pronounced scale separation of the dynamics. Specifically the following assumptions are introduced:

With these assumptions the plasma can be modelled in three parts: plasma transport, RF modulation and energy balance.

Plasma transport:

Radio Frequency Modulation:

Energy balance:

The model is completed by appropriate boundary conditions taking account of the sheath region. The first of these deals with ion velocity, where Bohm’s criterion is assumed:

. (35)

The sheath region is assumed to be depleted of electrons and unable to support drift or diffusion. Instead the current between bulk and electrodes is carried by the displacement current, with the sheath acting as a capacitor. With C_sh as the sheath capacitance and V as the electrode voltage this results in

. (36)

The electron depletion in the sheath also means that the plasma bulk is thermally insulated from the walls, i.e.

(37)

The above set of equations and boundary conditions define the EDDM.

Although the EDDM is derived from a fluid-dynamical model the method of implementation is very different and it is important that its behaviour should be as expected. Testing was carried out on a 1D parallel plate reactor containing an electro-positive (argon like) gas at a pressure of 250mTorr and with an applied RF voltage of 150V.

The electron and ion density for the fluid-dynamical model were compared with the plasma density for the EDDM (which assumes identical ion and electron densities). The results were qualitatively identical and quantitatively reasonable in the bulk region; of course, the EDDM makes no predictions in the sheath region.

Comparison of the averaged electron temperature showed that in each case the electron temperature increased towards the electrodes, in accordance with observation; this is due to the enhanced Ohmic heating in this region. In this case the quantitative agreement was very satisfactory.

The benefit of the EDDM is apparent from the fact that it required approximately two hundred times less computational time to converge.

The EDDM as incorporated in PHOENICS-CVD has been simplified to contain only a single ion species. This restriction simplifies the model and results in the solution of only four plasma equations in addition to the conventional transport equations. For most applications

there is a dominant carrier species which can reasonably be used to determine the parameters needed to define the plasma behaviour.

To facilitate the use of PHOENICS-CVD by CVD engineers with little, if any, CFD experience it has been equipped with a menu-driven user interface which makes problem set-up straightforward. The menu system is based on the standard PHOENICS menu system, with additional sections for the specification of options which are specific to CVD applications.

Special attention has been devoted to the requirements of the surface-to surface radiation model. This requires that all solid surfaces are covered by thermal zones, a task that is far from straightforward for the user. The menu system has therefore been designed to provide a starting point from which modifications and refinements can easily be made.

For more experienced users the alternative input using a PHOENICS Q1 file is still available, either as an alternative or a supplement to menu driven operation.

Additional data required by PHOENICS-CVD have already been mentioned in earlier sections; these are provided by the user in data files describing transport and thermodynamic properties, chemistry formulations and optical properties.

Extra features incorporated in PHOENICS-CVD are (i) the provision for the user to specify inflowing gas streams in conventional engineering units (sccm) rather than the standard PHOENICS mass flowrates, (ii) the introduction of special coding to allow the inclusion of shower plates as commonly used in CVD reactors, and (iii) the ability to model wafer batches approximately using a smaller number of wafers.

Throughout the ACCESS-CVD Project opportunities were taken to promote awareness of the work being done. This primarily took the form of conference and journal publications:

At the end of the Project a special two-day Seminar and Workshop was held in Wimbledon (September 28-29, 1995). During the Seminar presentations were made by the partners, describing all aspects of the Project: model development, validation and pilot implementation. The presentations were later published in a special edition of The PHOENICS Journal (December 1995, Volume 8, No 4), devoted exclusively to the ACCESS-CVD Project.