Numerical modelling of MPA-CVD reactors with the discontinuous Galerkin finite element method

In this article we develop a fully self consistent mathematical model describing the formation of a hydrogen plasma in a microwave power assisted chemical vapour deposition (MPA-CVD) reactor employed for the manufacture of synthetic diamond. The underlying multi-physics model includes constituent equations for the background gas mass average velocity, gas temperature, electromagnetic field energy and plasma density. The proposed mathematical model is numerically approximated based on exploiting the discontinuous Galerkin finite element method. We demonstrate the practical performance of this computational approach on a variety of CVD reactor geometries for a range of operating conditions.


Introduction
The remarkable optical, thermal, chemical and electronic properties of diamond, along with its extreme hardness and wear resistance, offer a material with great potential in a wide variety of industrial and scientific application areas. The synthesis of industrially grown diamond can be undertaken using two key approaches: high pressure high temperature (HPHT) and chemical vapour deposition (CVD). The former technique was originally developed in the mid 1950s, see [1], while CVD diamond synthesis was introduced in the mid 1980s, see [2], and the references cited therein. Since its introduction, the CVD diamond market has witnessed strong growth, due the strength, durability, stiffness and thermal conductivity properties of the resulting diamond product. Indeed, industrial CVD diamond synthesis offers the growth of high purity diamond, with low defect rates; this is crucial for the efficient manufacture of high quality products.
In this article we focus on the modelling of transport phenomena within a CVD reactor; for a detailed review of the key principles of effective CVD reactor design and its numerical simulation, we refer to [3]. Simulating CVD reactors represents an extremely challenging mathematical and computational task. Indeed, on the one hand, a critical feature is the design of the reactor geometry itself, in the sense that it must be able to support standing modes of the contained microwave field in such a manner that the location of peaks present in these standing modes must lie in the vicinity of the diamond substrate surface. These peaks lead to the ignition of a hydrogen plasma, which encourages the growth of diamond rather than graphite. Computational models which simulate the perturbation of the microwave field by the resulting plasma have been studied within the literature; for example, a finite difference technique was developed in [4]. This approach has subsequently been extended to include conservation of chemical species densities with application to the LIMHP reactor geometry, see [5]. An alternative numerical approach, based on employing the finite integration method, has been utilized in [6,7] for the simulation of the perturbed magnetic field applied to the AIXTRON reactor geometry, see below.
Secondly, the mathematical and computational challenge lies in modelling the combination of plasma physics, microwave field distribution, chemical reactions, and transport processes, as well as developing robust and efficient numerical discretization methods. We stress that the development of such models offer invaluable insight into the operating conditions during diamond growth, and moreover, allow for potential optim isation of the overall process, thereby, leading to improved reactor design for more effective diamond growth; work in this direction has been undertaken in [8], see, also, the references cited therein. A key mathematical issue is to balance the number of chemical species to be modelled, while ensuring an accurate description of the underlying CVD diamond synthesis process. Hassouni et al [9] present an excellent in depth review of hydrogen CVD plasma modelling; here, a model for the chemical species comprising hydrogen dissociation in CVD reactors is derived. This work highlights the complexity of plasma modelling due to the large number of chemical species involved in the many reactions. Moreover, this builds on prior research into CVD diamond synthesis modelling, and in particular, the coupling of the electric field with the plasma, as well as comparing full and simplified plasma physics models over a range of CVD operating conditions, see [5,10]. The complexity of more detailed methane plasma modelling is further discussed in [11]. However, we stress that the incorporation of complex flow phenomena in a given microwave power assisted (MPA-) CVD reactor model has not previously been taken into account; for example, [9] neglects the natural convection arising due to the Soret effect from strong thermal gradients, though it is stressed that these thermo-convective flow conditions may strongly affect the microwave field coupling configuration.
Stimulated by the work undertaken in [9], in this article we present a fully self consistent model for a hydrogen plasma, with simplified chemistry, within a MPA-CVD reactor. We point out that, as noted in [12,13], for example, the difference in the shape, position and optimal microwave power absorbed in a hydrogen plasma CVD reactor is not significantly different to the case when small quantities of methane are present. A key goal of this work is the incorporation of convective transport within the underlying model, together with the discretization of the proposed MPA-CVD model, whereby we simultaneously compute approximations to the microwave field, gas temperature, species composition, plasma density and background gas flow. Of particular interest, is the novel consideration of a background mass averaged gas flow; indeed, as noted above, previous work has simply considered diffusive transport phenomena only or modelled gas velocity by a velocity potential, see [14]. Here we compute the background gas mass average velocity field based on momentum conservation. In our numerical simulations, we compare computational results with those presented in the literature for designs inspired by the LIMHP and AIXTRON reactors. We then present a numerical experiment for the latter geometry which includes background gas flow; this demonstrates that the inclusion of convective transport can substantially alter the predictions of the MPA-CVD model. Throughout this article, the proposed MPA-CVD system of partial differential equations (PDEs) is discretized based on employing the discontinuous Galerkin (DG) finite element method (FEM). Broadly speaking DG schemes can be viewed as a hybrid between classical FEMs and finite volume methods (FVMs). Indeed, as in the FEM setting, the analytical solution is approximated based on employing polynomials of a given degree, defined over local element domains, but without the enforcement of any continuity constraints between neighbouring elements. Instead, elements are coupled based on employing numerical flux functions in an analogous manner to FVMs. Given this construction, DG schemes inherently possess a number of key advantages: improved stability for convection-dominated flows; local conservation property; simple treatment of meshes with hanging nodes when mesh adaptation is exploited; highly parallelisable due to limited interelement communication. For recent developments in the field, we refer to the recent monographs [15][16][17][18] and the references cited therein. In particular, DG methods are ideally suited to MPA-CVD applications, in the sense that they allow for the discretization of multi-physics problems in a unified manner. The definition of the underlying discretization is outlined in [19]; for brevity, we omit the details. Moreover, the implementation is based on employing the automatic code generation concepts outlined in [20] utilizing the AptoPy code developed in [19], together with the AptoFEM finite element package [21].
This article is structured as follows. In section 2 we introduce the self-consistent MPA-CVD reactor model. Section 3 then presents a summary of the proposed multi-physics system of PDEs, supplemented with appropriate boundary conditions. Numerical experiments highlighting the predictive capability of the proposed MPA-CVD model are presented in section 4; in particular, we consider CVD geometries inspired by the LIMHP and AIXTRON reactors. Finally in section 5 we summarise the work undertaken in this article and discuss potential future extensions.

Microwave power assisted chemical vapour deposition reactor model
The primary purpose of a CVD reactor model is to determine the shape, density and temperature of the ignited plasma above the diamond substrate surface. A self consistent description must couple the electric field propagation in the reactor cavity along with its perturbation resulting from the plasma. Furthermore, the transport phenomena arising from particle conservation of the low pressure hydrogen gas should also be taken into account. With this in mind, in this section we introduce the notion of conserving molecular and atomic hydrogen particle species' densities given their mass average flow field and heterogeneous chemistry due to chemical reaction. The time harmonic description of the CVD reactor's coupled electric field is then discussed for a medium whose electric permittivity is dependent on the characteristics of the reactor plasma. Finally, the electron particle conservation law is introduced, subject to the gain and loss of electrons due to ionisation and electron-ion recombination, respectively.

Mass averaged chemically reacting flow
In order to adequately model the position and fraction of dissociated hydrogen in a gas mixture of molecular hydrogen, the mass transport phenomena and heat of the gas mixture must be taken into account. In this section we consider flow for a multicomponent mixture comprising of n different chemically interacting species. The model presented accounts for conservation of each species' density, momentum and energy density. The theory of multicomponent flow is then applied to a binary gas mixture composed of molecular and atomic hydrogen. In this case, consideration will only be given to two components and the two chemical reactions enabling dissociation of molecular hydrogen. Here, the two components of atomic and molecular hydrogen are denoted by H and H 2 , respectively. The derivation of this model closely follows the theory of multicomponent flow outlined in [22].

Preliminaries of multicomponent flow
Consider the system of n interacting gaseous species; each gas species has density ρ i and velocity u i , = … i n 1, , . For each species, the rate of increase of mass must be balanced by the net rate of addition of mass to that volume by flow convection, molecular diffusion and chemical reaction. Thereby, the conservation of mass of component i gives where r i , = … i n 1, , , is the rate of mass production of component i due to chemical reaction. Summing over all components, = … i n 1, , , yields where we note that the total diffusive flux ∑ = = j 0 The conservation of component i in terms of molar concentration is given by where R i , = … i n 1, , , is the molar rate of mass production of component i due to chemical reaction, such that = r R M i i i . Summing equation (3) over all species and noting that moles are, in general, not conserved yields

Conservation of mass in a binary gas
Recalling the equations for mass and molar mass conservation, see (2) and (4), respectively, for the hydrogen binary gas, the diffusion fluxes j H and j H2 or * J H and * J H 2 must be expressed in terms of gradients of the concentration, temperature and pressure. Diffusion due to pressure gradients generally occurs only when the gradient in pressure is very large; typically, this is not the case in MPA-CVD reactors, and hence we neglect these terms. The diffusion flux will take into account concentration gradients of the mixture, as well as thermal terms. In general, for a multicomponent mixture, j i , = … i n 1, , , is determined by the Maxwell-Stefan equations; however, in the case of a binary gas mixture, the binary diffusion flux is given by Fick's law along with thermal terms. More specifically, for atomic hydrogen in a binary gas mixture of atomic and molecular hydrogen, we have that Here, HH 2 D is the binary diffusion coefficient for H in a mixture consisting of H and H 2 , D T H is the thermal diffusivity coefficient for H (note that = − D D T T H H 2 for a binary gas mixture of H and H 2 ), and T is the temperature of the mixture. The molar diffusion flux for a binary gas mixture of atomic and molecular hydrogen is given by the equivalent formulation of equation (5), namely, In essence the relationship between equations (5) and (6) is [19] for details. For the hydrogen binary gas, combining the equation of molar concentration conservation (4), with = i H, and the definition of the molar flux in equation (6) yields the conservation equation Substituting this into (7) yields the conservation of mass in terms of each species' molar concentrations and their gradients In order to define the density and concentration of a mixture, here the equations of state for an ideal gas are introduced: where P is the constant mean pressure of the system and R is the gas constant. Noting that ρ = Mc, ρ can be defined in terms of molar mass fractions comprising the mixture, i.e.
The binary gas mixture's molar mass fractions of atomic and molecular hydrogen can be expressed in terms of In the specific case of a binary gas mixture involving the dissociation of hydrogen, two simple reactions are assumed to take place, namely, The reaction rate for atomic hydrogen is given by , where R f is the forward molar rate of production by dissociation and R r is the reverse molar rate of production. Each of R f and R r are determined by the chemical kinetics of each reaction, depending on the concentration of each species H and H 2 along with the forward and reverse reaction rate coefficients, k f and k r , respectively. We refer to [22] regarding details of chemical kinetics in reaction terms; in the case of the rate of production of H by dissociation from H 2 we write where the factor of 2 accounts for 2H atoms produced per reaction. Here, we assume that the forward and reverse rate coefficients are equivalent for each reaction.

Conservation of momentum
Assuming the viscosity η of the mixture of n individual chemical species to be isotropic, the conservation of momentum is described by the Navier-Stokes equations, i.e.
Here, g is the acceleration due to gravity and the stress tensor τ is defined as where p is the pressure and I is the identity tensor. The viscosity of the mixture can be determined by the semi-empirical formula [23]; thereby, where η i is the viscosity of the ith component of the mixture, and Φ ij is given by In the case of the binary gas mixture of atomic and molecular hydrogen, we have that / , see [23].

Conservation of energy
The conservation of the energy density ρE of a multicomponent mixture is governed by where q is the energy flux of the mixture, Q is an energy source, ρ ⋅ u g is the rate of energy from external force acting on the gas mixture and E is the total energy of the system.
Here, E can also be expressed as = + E e u 1 2 2 , where e is the internal energy per unit mass, and u 1 2 2 is the kinetic energy per unit mass, where = ⋅ u u u 2 . We now consider the formulation of the mechanical energy equation: taking the product of the conservation of momentum equation (11) with u and rearranging with respect to the continuity equation, see [22], gives By subtracting the mechanical energy in equation (13) from the conservation of energy in equation (12) it can be shown that The conservation of energy can be expressed in terms of the system's temperature T and specific enthalpy h. The specific enthalpy for the multicomponent mixture can be expressed in terms of each of the specific enthalpies of the mixture's species, weighted by their respective mass fraction, The enthalpy of the system is composed of its internal energy in addition to the thermodynamic work done by the system on its adiabatic chamber, i.e. = + ρ h e p , where p is the pressure of the gas mixture. Substituting this expression for multicomponent enthalpy into the conservation of internal energy equation (14), it can be shown that where we have introduced the deviatoric component of τ, given In the CVD reactor model we assume the viscous dissipation effects in the CVD vacuum are negligible, since they are only important in flows with large velocity gradients; thereby, we set τ ∇ = u : 0 dev . Furthermore, given the ideal gas assumption of (8), i.e. that the gas flow in the CVD reactor is flowing in a system with constant mean pressure P, we set . Exploiting these assumptions, equation (15) reduces to An expression for the energy flux can be determined in terms of Fourier heat diffusion, where the Dufour effect of energy transport due to concentration gradients between chemical species is small and can be neglected [22]. Here, for multicomponent gas mixture temperature, The thermal conductivity κ can be expressed in terms of each chemical species thermal conductivity κ κ The term involving the spatial gradients of the species' specific enthalpies in the energy flux is often small and can be neglected, leaving the form of the conservation of energy in terms of enthalpy and temperature given by

Microwave field
The coupling of a microwave frequency electric field within the CVD reactor is key to the localised heating of the gas mixture above the diamond substrate. For a closed cavity whose walls act as perfect conductors, there are infinitely many harmonic electrical resonant modes at an equal number of corresponding frequencies which can be excited. Common commercially available magnetron devices operate at frequencies of 896 MHz and 2.45 GHz. A key aspect of CVD reactor design is to ensure that the reactor's geometry supports a given resonance mode at the frequency of the driving magnetron. Furthermore, it is favourable that the profile of the resonant electric field's amplitude be at its greatest magnitude above the substrate surface, ensuring efficient dissociation of hydrogen above where the diamond will be located. We model the electric field E, at angular frequency ω, in the MPA-CVD reactor plasma with the time harmonic formulation of Maxwell's equations, subject to Gauss's law for a dielectric medium, i.e. 1 2 where μ is the electric permeability and ε ε = − ′ σ ω j is the complex electric permittivity for conditivity σ and complex unit = − j 1. In the plasma ε ε = ′ p 0 , where ε 0 is the permittivity of free space, and the plasma conductivity is given by Here, ν me is the electron-neutral collision frequency and ω pe is the fundamental characteristic frequency of the plasma given by where q e is the charge of the electron, n e is the number density of electrons in the plasma, and m e is the rest mass of the electron.

Plasma
The source of electrons in the CVD reactor plasma can be modelled in a similar manner to the species of hydrogen, in the sense that their number density and internal energy should be conserved. As an alternative, we present a simplified model of electron momentum and energy conservation, along with electron generation by ionisation and loss by recombination with heavy ions according to the Maxwell-Boltzmann distribution. The following derivation summarises the plasma physics model presented in [24] and [25], see also [19]. The conservation of the number density of electrons in the MPA-CVD plasma requires consideration of at least seven species and their internal modes: H 2 , H, H + , H − , + H 2 , + H 3 , and e − . Under suitable simplifying assumptions of intermediate species, whose reactions are considered instantaneous, the complexity of modelling can be reduced to consider the five species: H 2 , H, H + , + H 3 , and − e , see [9]. Here, we have simplified further by modelling conservation of H and H 2 . We assume one ion species, H + , and thereby the number density of electrons is modelled by the ambipolar electron diffusion equation where n e is the number density of electrons, G is the gain of electrons by ionisation, L is the loss of electrons due to recombination and D a is the ambipolar diffusion coefficient.
Modelling the multitude of ionisation, excitation and recombination reactions of electrons with heavy species is an extremely challenging task. In this article, we consider the heuristic model employed in [26]; here, the collisions between the electrons and ions in the plasma, especially at low velocity, have a finite probability of resulting in their recombination into a neutral atom. This process is dependent on the collision cross section of electron recombination with the heavy ions Θ rec , the kinetic energy of the electron u e and the number density of the electrons and ions. This results in a loss of electrons in the plasma where Θ n u i rec e is the average collision frequency of the electrons with the ions. We can write (18) in terms of a recombination rate coefficient, which is a function of the electron temperature T e and the recombination energy E rec , i.e. The plasma is maintained in the steady state if the gain of electrons due to ionisation balances the losses due to diffusion and recombination. This process is dependent on the momentum exchange of electrons with the heavy species. In a similar fashion to the electrons lost by recombination, the electrons gained by ionisation is dependent on the average electron-heavy species ionisation cross section Θ i , average electron kinetic energy, the number density of electrons and the neutral heavy species n n , i.e. = Θ G n n u . e n i e This electron particle source term can also be written in terms of the rate coefficient k T i e ( ) and ionisation energy E i , namely, where, A i is the ionisation rate constant. A comprehensive model of the plasma chemistry occurring in the MPA-CVD reactor would account for mass and energy transport of the plasma electrons. To this end, the ambipolar diffusion equation (17) would be coupled to electron enthalpy conservation. This requires approximation of the electron thermal conductivity, as well as the rate coefficients of ionisation. The electrons would then be heated primarily by the applied microwave field; we refer again to [9] for a review of these models. In this work, we make the simplifying assumption that the electron temperature is strongly influenced by the time averaged magnitude of the applied electric field and the number density of the background gas of neutrals remains roughly constant such that for ionisation constant A 0 . We refer to [26] for appropriate choices of the constants R 0 and A 0 .

Ohmic heating
As a result of electron-neutral collisions arising from the electric field in a plasma, the time averaged collisional ohmic power absorbed by those electrons is defined by This source of heat is then applied in the conservation of energy of the multicomponent gas mixture heat source term Q of equation (12). Outside of the CVD reactor's vacuum region in which a plasma cannot be ignited, the standard Joule heating model is applied where σ is the electric conductivity of the medium.
For simplicity, here we assume that the energy absorbed in the discharge is entirely employed within the gas heating term, see (12); losses by radiation or energy carried away by excited particles to the walls, as a result of gas ionis ation and excitation of molecules, are neglected. However, we point out that in the case of rapid heating of nitrogen and air, the portion of energy transferred to gas heating typically lies within the range of 60%-90% of the total energy, depending on the operating conditions, see [27,28].

MPA-CVD reactor model problem
On the basis of the general mathematical model developed in the previous section, here we formally present the steadystate version of MPA-CVD reactor model problem posed in an azimuthally symmetric geometry. In particular, we write the model in a compact manner, which is appropriate for discretization purposes, as well as specifying appropriate boundary conditions for each variable in the underlying system of PDEs.
To this end, we write Ω to denote an open, azimuthally symmetric, bounded domain in 3 R with boundary ∂ Ω. An example of a typical CVD geometry is depicted in figure 1(a). Adopting a cylindrical coordinate system and assuming azimuthal symmetry, an axial slice at θ = 0 and r 0 ⩾ is taken yielding the bounded domain ⊂ Ω Ω * with boundary ∂ Ω * ; for simplicity of notation, we simply denote this two-dimensional slice by Ω. This domain is then subdivided into three subdomains characterising components of the CVD reactor such that Ω = Ω ∪ Ω ∪ Ω a q v . Here, Ω a is the subdomain filled with air at atmospheric pressure, Ω q is the fused silica window and Ω v is the vacuum chamber of the CVD reactor containing atomic and molecular hydrogen in which diamond growth occurs on the substrate surface.
The boundary of Ω is divided such that ∂ Ω = ∂ Ω surf ∪ ∂ Ω ∪ ∂ Ω ∪ ∂ Ω ∪ ∂ Ω ∪∂ Ω wall in out ant a xis . Here, ∂ Ω surf is the component of the boundary specifying the substrate surface on which the diamond is grown, ∂ Ω wall is the wall of the reactor, ∂ Ω in is the gas inlet, ∂ Ω out the gas outlet, ∂ Ω ant the microwave antenna which excites the electric field in the cavity and ∂ Ω axis the exterior boundary component which lies on the axis of symmetry located at = r 0. We write ∂ Ω ∂Ω , v q , and ∂ Ω a to denote the boundaries of Ω Ω , v q , and Ω a , respectively; furthermore, the interior subdomain interface boundaries are defined by Γ = ∂ Ω ∩ ∂ Ω aq a q and Γ = ∂ Ω ∩ ∂ Ω vq v q . On boundary Γ aq we define the unit normal vector pointing from Ω a to Ω q by n aq and from Ω q to Ω a by n qa . Similarly, on Γ vq we define the unit normal vector pointing from Ω v to Ω q by n vq and from Ω q to Ω v by n qv . The two dimensional slice of the CVD computational domain is shown in figure 1(b).

Conservation of momentum
Given an inlet gas flow profile u inlet , specifying a no slip boundary condition on the walls of the reactor and a Neumann condition on the symmetry axis, together with the outlet pipe, the multicomponent gas mixture momentum conservation equation takes the form:

Conservation of mass
Assuming that there is no presence of atomic hydrogen on the walls of the reactor and that the gas at the inlet is pure molecular hydrogen, the conservation law enforcing continuity of mass fraction is given by

Conservation of energy
The temperature within the reactor is modelled such that the walls and the inlet gas are held at room temperature T room and the diamond substrate surface is heated to T surface . The temperature and heat flux are required to be continuous across the interfaces between the vacuum, the quartz window and the air filled cavity. The energy of the CVD reactor is then model led by conserving temperature in the three different regions of Ω. Thereby, following [19], we consider the conservation law given by: The piecewise material parameters are given by where κ v , κ q and κ a are the thermal conductivities in the vacuum, quartz and air, respectively, and σ q and σ a are the electric conductivities of quartz and air, respectively. Here, κ v is given in terms of thermal conductivities of atomic κ H and molecular κ H2 hydrogen according to equation (16). Furthermore, we employ the thermal conductivity of fused silica κ q presented in [29].

Electron density
It is assumed that there are no free electrons present on all physical boundaries of the CVD reactor vacuum chamber and at the gas mixture inlet. Thereby, conservation of electron particle density is given by the ambipolar diffusion approximation

Electric field
The permeability of the gas mixture, quartz window and the air filled cavity are all assumed to be equivalent to the permeability of free space µ 0 . The complex permittivity ε ε = − ′ σ ω j is discontinuous across the subdomains of the reactor; namely,  [7,30,31]). The imposition of appropriate boundary conditions for modelling excitation of the electric field in the CVD reactor requires careful consideration; for a discussion of this topic for modelling antennae and waveguide feeds with finite element methods, we refer to [32]. In this work the tangential component of the incident electric field on the surface of the antenna is presumed to be known. To this end, enforcing perfect electric conductor boundary conditions on the reactor walls and excitation E ant from the antenna, the time harmonic for mulation of Maxwell's equations in the CVD reactor is given by Here, p is an additional (Lagrange multiplier) variable introduced in order that the divergence free condition of the electric field may be enforced when the underlying system of PDEs is discretized, see [33], for example.

Numerical experiments
In this section we present a series of numerical results to illustrate the predictive capability of the proposed MPA-CVD model within a variety of azimuthally symmetric CVD reactor geometries with different operating conditions. We stress that by restricting the geometry to a two-dimensional slice of an azimuthally symmetric CVD reactor, three dimensional affects, and in particular the modelling of non-azimuthally symmetric phenomena, cannot be taken into account; indeed, this forms part of our programme of future research. Here, the underlying system of nonlinear PDEs is numerically  approximated based on employing the DG finite element methodology; a key advantage of this approach is that it allows for the discretization of the underlying multi-physics problem in a unified manner. While, for reasons of brevity, the definition of the scheme are omitted (see [19] for details), we point out that the DG finite element formulation has been generated based on employing the symbolic algebra toolkit AptoPy [20]. AptoPy automatically generates the necessary Fortran code required by our inhouse software package AptoFEM [21] to assemble the underlying system of nonlinear equations, as well as the corre sponding Jacobi matrix for use within a (damped) Newton iteration. At each Newton step, the solution to the underlying linear system is computed using the MUltifrontal Massively Parallel (MUMPS) solver, see [34][35][36]. In order to compute the eigenpair solutions of the microwave cavity resonance problem, we employ the Arnoldi Package (ARPACK) [37]. For each test problem considered, the Triangle grid generation software package, see [38], has been exploited to generate the underlying finite element mesh. The chemical data relating to the dissociation of hydrogen and thermal properties of fused silica for all coefficients and parameters were obtained from the National Institute of Standards and Technology chemistry database [39,40].
In this section, we first present a numerical example to demonstrate results of the DG finite element approximation of the MPA-CVD reactor model in the simplified CVD geometry shown in figure 1. We then present results for reactor geometries inspired by the ellipsoidal AIXTRON reactor and the LIMHP reactor, both of which are summarised in [41]. Each reactor presented here is designed for operation with incident electromagnetic field operating at 2.45 GHz. For each geometry, the resonant electromagnetic field will be computed for the case of the empty, air filled, reactor. The peaks in the magnitude of the electric field solution provides an indication of the location of the ignited plasma when the hydrogen vacuum portion of the reactor is present. Given a satisfactory electromagnetic field configuration, we then proceed to compute a series of numerical solutions to the MPA-CVD model at hydrogen gas mixture pressures of 50 torr, 150 torr and 250 torr and input power of approximately 1 kW. The quantities of interest are the free electron density in the plasma n e , the number density of atomic hydrogen n H , the reactor temper ature T and the microwave power deposition in the plasma P ohm .

Example 1: a simple cylindrical MPA-CVD reactor
In our first example, we consider the geometry Ω depicted in figure 1(b). The reactor radius from the axis of symmetry to the outer wall is set at = r 12 cm reactor and the height of the chamber from the chamber floor to its ceiling is set at = h 24 cm reactor . The substrate surface has radius 3 cm and height 5 mm from the chamber floor. This configuration has been designed in such a manner to ensure that the transversal magnetic TM 022 harmonic mode is excited at a frequency of 2.45 GHz. In order to confine the ignited plasma to reside above the substrate surface, the base of the quartz window, whose width is set at 5 mm, is situated 6 cm from the floor of  the chamber. This position of the quartz window exploits the local maximum of the electric field TM 022 configuration at the base of the reactor. Finally, we note that the inlet gas nozzle is situated on the outer chamber wall 1 cm from the quartz window and has width 3 mm.
In order to numerically approximate the proposed MPA-CVD model, we construct the underlying finite element mesh such that there is a higher density of elements within the quartz window, near the gas inlet and outlet pipes, as well as in the vicinity of the substrate surface and microwave Figure 7. Example 2. The first generation LIMHP inspired reactor plasma shape operating at 1100 W antenna. The intention is to resolve the permit tivity transition between the vacuum and air filled cavity, the re-entrant corners of the geometry, and the regions of the domain where we expect a high density of electrons in the plasma to occur. The resulting mesh for this simple geometry consists of 44935 elements with 67040 interior and 725 boundary faces; see figure 2(a). Firstly, in figure 2(b) we show the magnitude of the electric field computed within the empty reactor; this yields a similar distribution to the empty cylindrical cavity resonator operating in a TM 022 mode, as expected. Thereby, we expect that the plasma will be ignited as intended above the diamond substrate surface. Furthermore, the lack of a peak in the electric field in the hydrogen gas mixture close to the quartz window should eliminate the possibility of damage to the window from unwanted heating.
On the basis of the mesh depicted in figure 2(a), the resulting DG finite element formulation for the MPA-CVD model in this configuration consists of 1358520 unknown degrees of freedom. This is split between: the velocity field variables u and p, 78225; the reactor temperature T, 134805; the species densities x H and n e , 33525 each; and the electric field variables E and p, 1078440. With this reactor design, we impose a zero inlet gas flow rate = u 0 inlet , such that diffusive transport effects dominate. In figure 3, we show the spatial distribution close to the substrate surface of the number density of atomic hydrogen n H , the number density of electrons n e , the temperature T, and the ohmic power P ohm , for power input 900 W and working pressures of 50 torr, 150 torr and 250 torr. The exploitation of the resonant electric field structure is clear with the shape and peak of the plasma density residing above the substrate surface. Combined with the ignited plasma's shape and position, this leads to the localised peak in the microwave power deposition above the substrate surface. At this power input the spatial gas temperature distribution is around the temperature required for efficient production of H from H 2 by dissociation, i.e. > T 3200 K. Figure 4 shows the temperature distribution in the quartz window, which does not exceed 650 K, well below its softening point of approximately 1343 K.
The electric field resonance property of the CVD reactor is affected by the permittivity of the generated plasma. Indeed, figure 5 demonstrates this variation in the electric permit tivity in the presence of the plasma. The attenuation of the electric field in the plasma is clear if we consider that it is approximately damped according to . In depth analysis of this coupling is undertaken in [42].

Example 2: the LIMHP Bell Jar reactor
In this section, the proposed MPA-CVD model is applied to the geometry depicted in figure 6(a) which approximates the first generation LIMHP bell jar reactor developed by LSPM Paris and Plassys analysed in [3]. In the geometry employed here, the reactor's 26 cm diameter Faraday cage houses a CVD chamber of height 35 cm containing a 10 cm diameter quartz  bell jar and 5 cm diameter substrate surface. The reactor is designed to operate at a mean pressure of 20 torr-120 torr, power of 0.6 kW-6 kW and microwave frequency of 2.45 GHz exciting the TM 023 cavity mode. The computational mesh generated for the LIMHP domain consists of 15836 elements with 23432 interior and 644 exterior faces. As before, regions of interest have a higher density of elements, such as the quartz bell jar, the microwave antenna and the substrate surface. In figure 6(b) we show the magnitude of the electric field in the empty air-filled LIMHP reactor. Here, we observe that there is a clear peak in the amplitude of the electric field near the substrate surface; however, a further peak is also present just above the top of the quartz bell jar. The primary peak in the electric field indicates that plasma ignition will occur as desired above the substrate surface, however, increasing the input power in the electric field does lead to a second plasma ball being generated at the top of the quartz bell jar. The first generation LIMHP reactor is known to exhibit this double plasma ball phenomenon [9].
Computational results for the MPA-CVD reactor model based on employing the DG finite element discretization are presented in figure 7 with very good agreement with those presented in [3]. In particular, we observe that the plasma density is distributed quite evenly across the substrate surface; moreover, its peak resides in close proximity to the carbon seed. However, as the power input and density increases, the formation of the second plasma ball at the top of quartz bell jar becomes prominent. The occurance of this discharge transition is explained in [42]. Due to both the peak in the electric field at the top of the quartz bell jar and the creeping temperature increase in the same region, a secondary plasma ball is generated. The detrimental effect of this secondary plasma ball is to shield the electric potential at the substrate surface, leading to suboptimal electron densities in the primary plasma. The attenuation experienced by the electric field is evident from the distribution of the dielectric constant shown in figure 8.

Example 3: the AIXTRON reactor
In this final example, we consider the performance of the proposed MPA-CVD model to a reactor design inspired by the ellipsoidal AIXTRON reactor depicted in figure 9(a). The ellipsoid design is intended to take advantage of its two focal points concentrating the microwave radiation at one point above the substrate surface. By ensuring resonance of the electric field at a frequency of 2.45 GHz with a TM 036 mode, the ellipse is constructed such that the minor radius = r 228 mm minor and major radius = r r major 4 3 minor . The base of the reactor is located 208 mm from its ellipsoid centre. The antenna is positioned in the cavity at the focal point . The diamond substrate surface has radius 2.5 cm and is raised to a height of 5 mm from the reactor floor. Finally, the gas inlet inner ring is located immediately next to the base of the substrate with inner radius = a 3.85 cm and outer radius = b 4 cm . The outlet outer ring has inner radius 5.8 cm and outer radius 6 cm. The computational mesh generated for this geometry consists of 34 243 triangular elements and 51 029 interior and 671 exterior faces; the resulting DG finite element approximation of MPA-CVD model consists of 1011 102 degrees of freedom.
The magnitude of the electric field computed in the empty air-filled cavity is presented in figure 9(b). Here, we observe that the peak in the resonant electric field mode lies on the diamond substrate surface; thereby, we expect the plasma predicted by our MPA-CVD model to reside as intended above the substrate. Furthermore, the small electric field amplitudes around the quartz bell jar should remove the risk of damage through unwanted microwave power deposition and hence heating.
Firstly, in this section we consider the case of a zero gas inlet flow rate, i.e. = u 0 inlet ; in this case, the spatial distributions of n H , n e , T, and P ohm are depicted figure 10.
Here, we observe that the shape of the plasma is improved by the AIXTRON inspired design over the simple geometry; indeed, the peaks in density of the plasma lie closer to the substrate surface and are more evenly spread. Furthermore, the temperature of the quartz bell jar does not exceed 400 K. We note that the results presented here for the AIXTRON inspired reactor are in agreement with those presented in [3].
Secondly, we consider the case when the inlet flow velocity is no longer zero. In this setting, we prescribe a standard parabolic coaxial Poiseuille inlet boundary condition at the inlet ring, namely, we set The scaling parameter u 0 allows for control of the inlet flow rate. The volume of gas pumped into the reactor is therefore In this case, we refine the computational mesh in the region near the inlet and outlet rings, see figure 11.  In figures 12 and 13 we present computational results for the DG finite element approximation of the proposed MPA-CVD model in the AIXTRON inspired reactor with inlet gas flow rate of 10 l − min 1 at working pressure 250 torr and power input of 1100 W. In particular, in figure 12 we overlay the temperature field with the background gas flow streamlines above the substrate surface. Here, the width of the streamlines is weighted by the normalised magnitude of the gas velocity. Furthermore, figure 13 shows the atomic hydrogen density, electron density, temperature and ohmic power fields above the substrate.
We can determine from these results that the gas flow rate affects both the temperature and atomic hydrogen density. Indeed, the convective gas transport leads to a sharper peak in the gas temperature and atomic hydrogen density above the substrate surface. The spatial distribution of the electron density appears to be almost entirely unaffected by convective transport. The resulting increase in the plasma density is due to the change in gas temperature and hence the electron neutral collision frequency ν me and ambipolar diffusion coefficient D a .

Conclusions and future work
In this article we have proposed a fully self consistent mathematical model for approximating the physical phenomena which arise in a MPA-CVD reactor for application in diamond manufacturing. In particular, we have considered the case of a hydrogen gas plasma discharge; the resulting system of nonlinear PDEs underpinning this model represents an extremely challenging computational task. Here, we have exploited the DG finite element method, based on employing the AptoPy and AptoFEM software packages, see [21] and [19], respectively. Computational examples have been presented which clearly highlight the predictive capability of the proposed approach; indeed, simulations for MPA-CVD reactor designs inspired by the LIMHP reactor and the AIXTRON reactor demonstrate excellent agreement with those presented in [3]. Moreover, the effect of varying the gas inlet velocity has also been studied for this latter reactor geometry.
Future work will be devoted to improving the predictive capability of the underlying plasma physics model. For example, on the basis of the work undertaken in [9], the binary gas phase model discussed in this article can be generalised to a full multiphase model which includes, in its simplest form, the species of H, H 2 , and e. In the case of more comprehensive models of hydrogen plasma chemistry, we need to consider more than a single ion species; thereby, the ambipolar diffusion approximation would not be sufficient. This requires careful consideration of the large number of reactions between the numerous species of hydrogen dissociation. Considering the electrons in the plasma discharge as a separate chemical species allows for a more accurate description of their convective and diffusive mass and energy transport. However, this requires a more precise representation of their generation and loss through ionisation, excitation and recombination. In addition, improvements to the specification of the boundary condition modelling the excitation of the electric field will be studied. We also plan to consider the extension to three-dimensional CVD reactor geometries, in order to model non-azimuthally symmetric phenomena.