## 3 Formulation and Numerical Modeling
### Governing Equations  
In order to consider the microwave heating of plasmas in the MPJCVD system, the microwave plasma interface is employed as implemented in COMSOL. The electron dy  
Figure 3: Research flow for the optimization of the MPJCVD chamber configuration with the FEM modeling for diamond film fabrication and characterization. namics are described by a pair of DD equations, one equation, for the electron density and one for the electron energy density, as follows [43]:  
\[\frac{\partial n_{e}}{\partial t}+\
abla\cdot\Gamma_{e}=R_{e}, \tag{1}\]  
\[\frac{\partial n_{e}}{\partial t}+\
abla\cdot\Gamma_{e}+\mathbf{E}\cdot\Gamma_{ e}=R_{e}+\frac{Q_{rh}}{|q|}, \tag{2}\]  
Here, \(n_{e}\) is the electron density, \(n_{e}\) is electron energy density, \(\Gamma_{e}\) is the electron flux, \(\Gamma_{e}\) is the electron energy flux, \(\mathbf{E}\) is the quasi-static electric field, \(R_{e}\) represents the electron source, \(R_{e}\) is the energy loss due to inelastic collisions, \(Q_{rh}\) is the power transfer from the electromagnetic field to the electrons, and \(q\) is the electron charge. The electron flux \(\Gamma_{e}\) and the electron energy flux \(\Gamma_{e}\) are defined explicitly as  
\[\Gamma_{e}=[-n_{e}(\mu_{e}\cdot\mathbf{E})-D_{e}\cdot\
abla n_{e}] \tag{3}\]  
\[\Gamma_{e}=[-n_{e}(\mu_{e}\cdot\mathbf{E})-D_{e}\cdot\
abla n_{e}] \tag{4}\]  
where \(\mu_{e}\) represents the electron mobility, \(\mu_{e}\) is the electron energy mobility, \(D_{e}\) is the electron diffusivity, and \(D_{e}\) is the electron energy diffusivity. They are related and calculated by  
\[D_{e}=\mu_{e}T_{e}, \tag{5}\]  
\[\mu_{e}=\frac{5}{3}\mu_{e}, \tag{6}\]  
\[D_{e}=\mu_{e}T_{e}. \tag{7}\]  
where \(T_{e}\) is the electron temperature. The source coefficients \(R_{e}\) and \(R_{e}\) are determined by the plasma chemistry and are written using rate coefficients. Suppose that there are \(M\) reactions that contribute to the growth or decay of electron density and \(P\) inelastic electron-neutral collisions. The electron source \(R_{e}\) and the energy loss due to inelastic collisions \(R_{e}\) are calculated respectively by  
\[R_{e}=\sum_{i=1}^{M}x_{i}k_{i}N_{n}n_{e}, \tag{8}\]  
\[R_{e}=\sum_{i=1}^{P}x_{i}k_{i}N_{n}n_{e}\Delta e_{i}, \tag{9}\]  
where \(i\) is the reaction index, \(x_{i}\) is the mole fraction of the target species for reaction \(i\), \(k_{i}\) is the rate coefficient for reaction \(i\), \(N_{n}\) is the total neutral number density, and \(\Delta e_{i}\) is the energy loss due to reaction \(i\). The rate coefficients can be computed from cross section data by the following integral:  
\[\frac{\partial n_{e}}{\partial t}+\
abla\cdot\Gamma_{e}=R_{e}, \tag{1}\]0  
where \(\epsilon\) is the energy, \(m_{e}\) is the electron mass, \(\sigma_{i}\) is the collision cross section for reaction \(i\),\(f\) is the EEDF, which we assume is Maxwellian as mentioned above. The electron source and inelastic energy loss are automatically computed by the COMSOL multiphysics interface once the chemical species and reactions are chosen. For non-electron species, the mass fraction \(w_{k}\) of each species \(k\) can be solved by the convection-diffusion equation:  
\[\frac{\partial n_{e}}{\partial t}+\
abla\cdot\Gamma_{e}=R_{e}, \tag{1}\]1The electrostatic field is calculated by solving the Poisson equation for the electrostatic or scalar potential V and taking a minus gradient of it, i.e., \(-\
abla\)V:  
\[\frac{\partial n_{e}}{\partial t}+\
abla\cdot\Gamma_{e}=R_{e}, \tag{1}\]2  
where \(\varepsilon_{0}\) is the permittivity of vacuum and \(\varepsilon_{r}\) is the relative permittivity of materials in the simulation domain, which usually has a spatial dependence. The space charge density \(\rho\) is computed based on the plasma chemistry using the equation below:  
\[\frac{\partial n_{e}}{\partial t}+\
abla\cdot\Gamma_{e}=R_{e}, \tag{1}\]3  
where \(Z_{k}\) is the charge and \(n_{k}\) is the number density for species \(k\).