Phenomenological Single-Particle
Modeling of Reactive Transport
in Semiconductor Processing

3.5 Extended Knudsen Diffusion

The standard Knudsen diffusion calculation discussed in Section 3.2 requires substantial approximations. Namely, the sticking coefficients are assumed to be low, and the features are assumed to be long, such that the diffusion process can be approximated by the diffusivity of infinite cylinders or trenches. Therefore, local contributions to the flux due to visibility of the source (c.f. Section 2.3.2) are neglected. Although Clausing has calculated corrections in transmission probabilities for short cylinders [114, 115], the issue of direct source flux contributions remains unaddressed.

The hitherto discussed formulation of Knudsen diffusivity including the view factors enables an extended calculation overcoming these limitations for short cylinders and including the direct flux in a straightforward manner. In essence, the mass balance procedure follows the same derivation that Eq. (3.15), however, the preferential transport direction $z$ is no longer assumed to extend from $-\mathrm{\infty }$ to $\mathrm{\infty }$. Instead, the same mass balance integrals to calculate ${\mathrm{\Gamma }}_{\mathrm{cross}}$ from Eq. (3.12) are now restricted to the range $0$ to $L$, as illustrated in Fig. 3.7. This restriction now enables considering the contributions due to the direct flux from the source area $A_{\mathrm{source}}$.

In contrast to the bottom-up approach discussed in Section 2.3.2, the direct flux is not computed at each surface element. Instead, it is accumulated at the entire cross-section through the $N^{-}$ element of the net flux balance in Eq. (3.12). It now reads

$$\begin{array}{}\text{(3.45)}& N^{-}={\mathrm{\Gamma }}_{\mathrm{source}}{A}_{\mathrm{source}}{F}_{\mathrm{source}-\mathrm{cross}}(0-z)+{\int }_0^z(1-\beta (z^{\prime })){\mathrm{\Gamma }}_{\mathrm{imp}}(z^{\prime })F_{d{A}^{\prime }-A}(z^{\prime }-z)d{A}^{\prime }.\end{array}$$

In Eq. (3.45), $F_{\mathrm{source}-\mathrm{cross}}$ is the finite-finite view factor between $A_{\mathrm{source}}$ and $A_{\mathrm{cross}}$ [98]. This element is not multiplied by a sticking coefficient since it is assumed that the source is fully emitting (i.e., ${\beta }_{\mathrm{source}}=0$). Similarly, the $N^{+}$ term must consider the flux due to reflections from the bottom area $A_{\mathrm{bottom}}$:

$$\begin{array}{}\text{(3.46)}& N^{+}=(1-\beta (L)){\mathrm{\Gamma }}_{\mathrm{imp}}(L)A_{\mathrm{bottom}}{F}_{\mathrm{bottom}-\mathrm{cross}}(L-z)+{\int }_z^L(1-\beta (z^{\prime })){\mathrm{\Gamma }}_{\mathrm{imp}}(z^{\prime })F_{d{A}^{\prime }-A}(z^{\prime }-z)d{A}^{\prime }\end{array}$$

Just as in Section 3.2, in order to avoid an integral equation, a Taylor expansion of the concentration similar to Eq. (3.16) is required. There are, however, subtle differences. The lack of infinite dimensions means that even-order terms do not cancel from symmetry considerations. Instead, all terms above the first derivative are directly truncated, and the term involving $n(z)$ cannot be disregarded.

For simplicity, similar considerations to those made in Eq. (3.31) are made. First, a constant value of $\beta$ is assumed. Also, all terms are normalized with respect to ${\mathrm{\Gamma }}_{\mathrm{source}}$, therefore, the Taylor expansion can be performed over ${\hat{\mathrm{\Gamma }}}_{\mathrm{imp}}(z)$ instead of $n$, as they are equivalent after normalization (c.f. Eq. (3.7)). Additionally, in a similar vein to Eq. (2.7), a normalized cross-sectional flux is defined as:

$$\begin{array}{}\text{(3.47)}& {\hat{\mathrm{\Gamma }}}_{\mathrm{cross}}(z)=\frac{{\mathrm{\Gamma }}_{\mathrm{cross}}(z)}{{\mathrm{\Gamma }}_{\mathrm{source}}}\end{array}$$

With all these considerations, Eqs. (3.45) and (3.46) are combined with the net flux from Eqs. (3.12) and (3.47) to obtain

$$\begin{array}{}\text{(3.48)}& \begin{array}{c}{\hat{\mathrm{\Gamma }}}_{\mathrm{cross}}(z)=\\ \frac{A_{\mathrm{source}}}{A_{\mathrm{cross}}(z)}{F}_{\mathrm{source}-\mathrm{cross}}(0-z)-(1-\beta )\frac{A_{\mathrm{bottom}}}{A_{\mathrm{cross}}(z)}{F}_{\mathrm{bottom}-\mathrm{cross}}(L-z){\hat{\mathrm{\Gamma }}}_{\mathrm{imp}}(z)\\ +(1-\beta )\frac{1}{A_{\mathrm{cross}}(z)}[{\int }_0^z{F}_{d{A}^{\prime }-A}(z^{\prime }-z)-{\int }_z^L{F}_{d{A}^{\prime }-A}(z^{\prime }-z)d{A}^{\prime }]{\hat{\mathrm{\Gamma }}}_{\mathrm{imp}}(z)\\ +(1-\beta )\frac{1}{A_{\mathrm{cross}}(z)}[{\int }_0^z(z^{\prime }-z)F_{d{A}^{\prime }-A}(z^{\prime }-z)-{\int }_z^L(z^{\prime }-z)F_{d{A}^{\prime }-A}(z^{\prime }-z)d{A}^{\prime }]\frac{d{\hat{\mathrm{\Gamma }}}_{\mathrm{imp}}(z)}{dz}\\ -(1-\beta )\frac{A_{\mathrm{bottom}}}{A_{\mathrm{cross}}(z)}(L-z)F_{\mathrm{bottom}-\mathrm{cross}}(L-z)\frac{d{\hat{\mathrm{\Gamma }}}_{\mathrm{imp}}(z)}{dz}.\end{array}\end{array}$$ Using the same conventions, the steady-state conservation of mass can be written as:

$$\begin{array}{}\text{(3.49)}& \frac{d{\hat{\mathrm{\Gamma }}}_{\mathrm{cross}}}{dz}=-\bar{s}\beta {\hat{\mathrm{\Gamma }}}_{\mathrm{imp}}(z)\end{array}$$ The BCs are then equivalent to those from Eqs. (3.29) and (3.30):

$$\begin{array}{}\text{(3.50)}& {\hat{\mathrm{\Gamma }}}_{\mathrm{imp}}(0)=& 1\text{(3.51)}& {\hat{\mathrm{\Gamma }}}_{\mathrm{cross}}(L)=& \beta {\hat{\mathrm{\Gamma }}}_{\mathrm{imp}}(L)\end{array}$$

It is important to note that, instead of a second-order ODE in $n$, the extended Knudsen diffusion is a system of coupled first-order ODEs for ${\hat{\mathrm{\Gamma }}}_{\mathrm{imp}}$ and ${\hat{\mathrm{\Gamma }}}_{\mathrm{cross}}$. Although it is indeed not an integral equation, it still requires the pre-computation of several integrals which might have fairly complex forms. This, in combination with the complexity of Eq. (3.48), makes the use of numerical methods necessary.

To investigate the consequences of this extended calculation, it is evaluated for a finite cylinder of constant diameter. The same differential-finite view factor from Eq. (3.23) is used, as well as the disk to parallel coaxial disk view factor, necessary for the terms involving both $A_{\mathrm{source}}$ and $A_{\mathrm{bottom}}$ [107, 128]

$$\begin{array}{}\text{(3.52)}& F_{\mathrm{disk}-\mathrm{disk}}=\frac{1}{2}[R-\sqrt{R^2-4}],\end{array}$$ where

$$\begin{array}{}\text{(3.53)}& R=\frac{2\frac{d^2}{{z^{\prime }}^2}+4}{\frac{d^2}{{z^{\prime }}^2}},\end{array}$$ such that both disks are separated by $z^{\prime }$ and have equal diameters $d$. The resulting ODEs after computing the involved integrals in Eq. (3.48) have a closed-form expression, however, they are omitted for brevity. The numerical solution of the ODE system is computed using Mathematica [126]. It is compared to both the solution of the standard Knudsen diffusion from Eq. (3.31) and that obtained with a radiosity framework [108], shown in Fig. 3.8.

In Fig. 3.8, it can be seen that the extended Knudsen diffusion follows very closely the curve obtained using the radiosity framework. Since the latter evaluates the integral equation and has been validated with a Monte Carlo simulation [108], it can be considered the exact result for a cylinder. For cylinders with low $\beta$, all flux calculations yield very similar results which is expected from the characteristics of Knudsen diffusion discussed so far. Interestingly, qualitatively similar results among all calculations are also obtained for the lowest AR cylinder, even though the standard Knudsen diffusion has been constructed for long cylinders. This is evidence that, for relatively short cylinders, the important physical phenomena are captured by the boundary conditions instead of the diffusivity.

In the situation of high AR and high $\beta$, standard Knudsen diffusion deviates more notably from both extended Knudsen diffusion and the radiosity framework. This is due to the lack of direct flux in standard Knudsen diffusion, thus the exponential decline of Eq. (3.31) is the only dominating factor. Naturally, in situations of higher $\beta$, attention must be placed on the role of the direct flux. A methodology to partially recover the effects of the direct flux while still using the standard Knudsen diffusion is discussed in Section 3.6.2.

Nonetheless, the success of the extended Knudsen diffusion should be carefully interpreted. The qualitative agreement shown in Fig. 3.8 is an indication of the small relative error between the extended diffusive calculations and the exact radiosity framework. However, since ${\hat{\mathrm{\Gamma }}}_{\mathrm{imp}}$ is already normalized to ${\mathrm{\Gamma }}_{\mathrm{source}}$, it is in fact more useful to focus on the absolute error in ${\hat{\mathrm{\Gamma }}}_{\mathrm{imp}}$ to obtain the deviation in units of ${\mathrm{\Gamma }}_{\mathrm{source}}$. This is shown in Fig. 3.9 for the case of a cylinder of $d=1\text{a.u.}$ and $\text{AR}=100$. There, it can be seen that the error is more pronounced at the top of the cylinder for higher values of $\beta$. Such error is likely a consequence of the boundary condition in Eqs. (3.29) and (3.50) imposing maximum flux at the top. For higher values of $\beta$, the reduction in flux due to the $90{}^{\circ }$ inclination of the cylinder wall with respect to the source plane is an important factor which is not captured.

Additionally, the difference in absolute error between the standard and extended Knudsen diffusion is minimal. This is more apparent when one of the results from Fig. 3.8 is evaluated without the logarithmic scale. This is shown in Fig. 3.10 for the cylinder with $\text{AR}=100$ and $\beta =1\mathrm{\%}$. In this linear plot, it is clear that all curves are very similar qualitatively. Therefore, the standard Knudsen approach is more than adequate for many applications, as long as the phenomenological parameters are sufficiently adjusted.