A set of numerical tools for the analysis and dynamic dimension reduction of chemical vapor and atomic layer deposition (ALD) surface reaction models is developed in this work. The approach is based on a two-step process where in the first, the chemical species surface balance dynamic equations are factored to effectively decouple the (nonlinear) reaction rates, a process that eliminates redundant dynamic modes and that identifies conserved quantities. If successful, the second phase is implemented to factor out redundant dynamic modes when species relatively minor in concentration are omitted; if unsuccessful, the technique points to potential model structural problems. An alumina ALD process is used for an example consisting of 19 reactions and 23 surface and gas-phase species. Using the approach developed, the model is reduced by nineteen modes to a four-dimensional dynamic system without any knowledge of the reaction rate values. Results are interpreted in the context of potential model validation studies.

## I. INTRODUCTION

There is a long history associated with the mathematical analysis of chemical reaction dynamics to generate reduced-order models, particularly in homogeneous reaction systems.^{1} This is particularly true of combustion reaction dynamics research, where the scale and complexity of the reaction networks (RN) has motivated a number of important studies.^{2–4} By comparison, heterogeneous reacting systems have received less attention, especially in the realm of thin-film processing.

Effective modeling of heterogeneous reaction dynamics is crucial to atomic layer deposition (ALD) high-throughput system process design and optimization.^{5,6} In ALD, thin films are deposited in discrete (usually submonolayer) increments by exposing the growth surface to the gas-phase reactive precursors in an alternating manner. Unlike chemical vapor deposition, ALD does not possess steady state growth modes; instead, continuous ALD operation consists of limit-cycle behavior.

Generating accurate experimental data useful for studying intrinsic ALD reaction rates is challenging because of the dynamic nature of the process, the wide time-scale range of the dynamic processes, and the small surface concentrations of some of the intermediate reaction species. Therefore, researchers have turned to quantum chemical computations to model and compare ALD reaction mechanisms^{7–10} and to validate those predictions against the data that are available. Because these methods produce static information regarding reaction energetics and transition state configurations, the transitions between surface states can only be described in terms of equilibrium relations. However, conventional transition-state theory (CTST)^{11} does provide rate expressions for the activated surface reactions, and rate equations for precursor species adsorption and desorption processes likewise can be formulated.

Given the modeling situation described, we should expect the overall ALD RN model to consist of a differential-algebraic equation (DAE) system; a key objective of this work is to develop rational methods to formulate ALD reaction kinetics models in the form of a well-posed DAE system. By writing all reactions in terms of their net-forward rates, including equilibrium processes by adding a time constant *ϵ* ≪ 1 s, the pure differential-equation system then can be factored to decouple the reaction terms.^{1,12} Success of this factorization procedure indicates that the outer solution^{13} to this singular perturbation problem (found for *ϵ* → 0) retains the dynamics required to accurately model the deposition process. A number of reaction-independent modes also are normally produced by this process, reflecting the conserved quantities of the process and the elimination of redundant dynamic modes. A further reduction of dynamic dimension then can be gained by repeating the factorization process for the species balances associated with the slow and conserved modes, focusing only of species identified as being major relative to (minor) intermediate surface species. This phase of the reduction process thus identifies combinations of finite-rate processes that approximate new equilibrium relationships. The net result is a procedure where a minimal-order dynamic ALD reaction model is generated, a computational process during which potential model structural deficiencies are identified.

## II. PROTOTYPE DEPOSITION SYSTEM

To describe the model reduction and analysis process we developed, consider the simplified gas/surface RN and the net-forward rates *f _{i}* associated with each reaction shown in Fig. 1, where M and D are gas-phase monomer and dimer species with number concentrations [

*M*] and [

*D*] in m

^{−3}(e.g., for alumina ALD, M represents the trimethylaluminum monomer and D its gas-phase dimer). A is an adsorbed surface species with number concentration [

*A*] m

^{−2}while B represents the element of precursor M that is incorporated into the bulk film and A

^{‡}the transition state of the final irreversible reaction. Note that any by-products of the irreversible deposition reaction are omitted from this example. Sites open to adsorption S are consumed during adsorption of M and regenerated as bulk film B is created. The concentration of B [

*B*] also has a unit m

^{−2}; its value can grow infinitely large because [

*B*] represents the total number of atoms deposited per unit surface area.

The first three reactions of Fig. 1 are reversible, and the first of these is modeled as being barrierless in the forward direction, and so no transition state is defined for the adsorption process. Reaction *f*_{0} is a gas phase reaction while the remainder is considered (in terms of their dimensions) surface reactions.

For this RN with *σ* = surface area/volume, a molar balance on each of the six species gives

subject to the specified initial condition

Conceptually, we can think of this as an isolated system of volume *V* and reaction surface area *σV*. For simplicity, we take *σ* = 1 m^{−1} and generally omit references to units in the following discussion of this prototype system.

We define the forward rates of dimerization, monomer adsorption, and subsequent surface reaction processes as

where *ϵ* is the time constant (this time constant has units of s and may be an artificial construct, as will be described later) corresponding the reactions (3) and (5) dynamically relaxing to chemical equilibrium defined by *g*_{0} = 0 and *g*_{2} = 0. It is possible to solve Eqs. (1)–(6) directly for finite *ϵ* using a numerical integration technique suitable for stiff problems, although this solution can only be considered the true solution when *ϵ* = 0 (but which is impossible to compute directly). However, we first examine the role reaction factorization^{12} has in removing redundant dynamic modes before performing the numerical integration for *ϵ* > 0.

A systematic approach to developing a dynamic model of the surface reaction processes is based on performing a Gauss–Jordan factorization^{12,16} on stoichiometry matrix **S** defined in Eq. (1) to find

or

where the *x _{i}* are the current dynamic states and

*w*are constants. The factorization procedure is based on finding the correct linear combinations of each species balance to effectively decouple the reaction rates in Eq. (7), e.g., add twice the first balance (for D) to the second (for M) to eliminate the influence of

_{j}*f*

_{0}from the second differential equation.

The results of this relatively simple operation are far-reaching. First, we observe that the conserved (instantaneous) quantities *w*_{0} and *w*_{1},

are constants representing in Eq. (8) the total number density of adsorption sites, and in Eq. (9), the total number of deposition element atoms in this isolated system. Note that the constant *σ* = 1 m^{−1} premultiplies each surface concentration term to maintain consistent dimensions in Eq. (9). Given the clear physical significance of *w*_{0} and *w*_{1}, we conclude that *the factorization procedure can reveal essential features of the reaction chemistry.*

Computational results for *ϵ* = 0.1 and *K*_{0} = *K*_{1} = *K*_{2} = 1, *k*_{1} = 2, and *k*_{3} = 1 are presented in the left plot of Fig. 2; note that solutions to Eqs. (1) and (7) are identical (and must be). For this simulation, the specified initial condition is [*M*]_{0} = [*S*]_{0} = 1 and all remaining species concentrations are zero. We observe in Fig. 2(left) that the monomer M concentration drops immediately, and in response, the dimer D concentration rises. After this extremely quick dynamic rearrangement, the monomer M and dimer D begin to be consumed by the deposition reactions, resulting in a quick rise in the adsorbed species concentration [*A*], followed by a slower tailing off of [*A*]. The concentration of deposited species [*B*] rises throughout the simulation period from its initial state of zero; the concentration of open deposition sites [*S*] drops quickly at the outset as [*A*] and [*A*^{‡}] (the latter not shown) increase, and then rises, ultimately to return to its initial value as the precursor species are consumed.

### A. Singular perturbation solution

If the dimerization reaction is considered fast relative to the adsorption and solidification reactions and because of the requirements of CTST,^{11} the true dynamics are described by a differential-algebraic system of equations that result in the limit of *ϵ* → 0. Setting *ϵ* = 0 in (1, 3–6) directly gives *g*_{0} = *g*_{2} = 0 and one differential equation each for $[S\u0307]$ and $[B\u0307]$—an underdetermined system for the six species of this RN. This illustrates the challenge of writing a well-posed set of DAE without the factorization procedure we pose, an issue relevant to a wide range of CVD and ALD systems.

Proceeding from the factored system (7), however, the outer (slow time) solution of this singular perturbation problem for *ϵ* → 0 is found by numerically integrating the semiexplicit DAE system

subject to **c**^{0} consistent with $x10,\u2009x30$ and the algebraic equations (12)–(15). A solution corresponding to the specified initial conditions (2) is shown in Fig. 2(center), where the effect of projecting Eq. (2) onto Eqs. (12)–(15) can be seen as the instantaneous jump of [*D*] to a nonzero value and of [*M*] to a value less than that of Eq. (2).

Setting *ϵ* = 0 eliminates the rapid equilibration phase of the dimerization reaction; after that point, the dynamics are essentially equivalent to the previous system featuring four dynamic modes. An alternative view of the jump from the specified initial condition to an initial state consistent with Eqs. (12)–(15) is seen as the red line segment of Fig. 3 connecting **c**_{o} to **c**^{0}—conceptually, both of these points correspond to *t* = 0. For *t* > 0, the trajectory follows the red squares, which lie on the two-dimensional manifold $Q$ in the three-dimensional species space shown in this diagram; this manifold is defined by values ([*M*], [*D*], and [*A*]) such that Eqs. (12)–(15) are satisfied.

### B. Further dynamic dimension reduction

For many deposition systems, the concentration of transition-state species will be vanishingly small and so $[A]\u226b[A\u2021]$; furthermore, it is common to find that the adsorbed species occupies a small fraction of the potential adsorption sites and so $[M],[D]\u226b\sigma [A]$. In this situation

A physical interpretation of Eq. (16) is that it defines the pseudoequilibrium between the adsorption and reaction steps of a typical deposition process and that relationships of this type are used to define the Arrhenius plots characterizing mass-transfer versus kinetically limited CVD regimes. A detailed discussion of these operating regimes in the context of dynamic ALD processes is given in Adomaitis.^{14}

At this point, the dynamical behavior of the originally six-dimensional system (1) has been reduced to one in the new coordinate; the single differential equation in time plus the five algebraic equations constitute a well-posed problem because it is relatively easy to combine Eq. (10) with Eqs. (12)–(14) and (16) to ultimately find

where the other species concentrations are then computed from Eqs. (12) to (15) and (16). It is easy to see in Eq. (17) that [*M*] → 0 as *t* → *∞* and that the steady-state solution [*M*] = 0 is unique and asymptotically stable.

The solution to this reduced-order system corresponding to specified initial condition (2) is shown in Fig. 2, right. As with the two-dimensional system, [*D*] and [*M*] jump immediately to their pseudoequilibrium values at *t* = 0. Additionally, we observe that the initial adsorption dynamics of [*A*] and [*S*] are no longer present; they too immediately take on their pseudoequilibrium values. Most importantly, however, is that the deposition rate of B is essentially unaffected by the dynamic reduction process. This fact is vitally important for accurate modeling of the deposition process dynamics.

As with the other cases, we also represent the fully reduced model in Fig. 3; the jump from the specified initial condition to an initial state consistent with Eqs. (12)–(15) and *f*_{1} − *f*_{3} = 0 is seen as the green line segment connecting **c**_{o} to **c**^{0}. After this initial jump, the trajectory follows the green, one-dimensional manifold $QM$ where the points in time are shown as the yellow circles. The pseudoequilibrium manifold $QM$ for this case is defined by

To conclude the analysis of this prototype deposition RN, we consider one final view of the multiple time-scale features of Fig. 3: following the dotted black curve corresponding a finite-*ϵ* solution from the specified initial state **c**_{o}, we observe the rapid relaxation from the initial state (2) to the pseudoequilibrium condition $Q$ of the gas-phase dimerization reaction (a fixed time-step size was used and so closer spacing of points indicates slower dynamics). The dynamics of adsorption of species M to form A govern the next phase of this transient process. After the adsorption and irreversible deposition (solidification) process forming B reach a pseudoequilibrium, the deposition process enters its final dynamic phase where the gas-phase reactants are more slowly depleted, resulting in a final concentration of deposited material [*B*]^{∞} with all other species concentrations approaching zero along manifold $QM$. Overall, the multiple time-scales of this problem manifest themselves as a contraction of the initial dynamics in the three-dimensional ([*M*], [*D*], and [*A*])-species concentration space to a two-dimensional equilibrium manifold $Q$, then to the one-dimensional manifold $QM$, to the zero-dimensional equilibrium **c**^{∞}. The final, zero-dimensional state corresponds to steady CVD operation or an infinitely long ALD exposure.

### C. Systematic factorization procedure

While hand-calculations are straightforward to implement for the reaction factorization/dynamic mode reduction process on our prototype RN, a more systematic procedure is needed for RN of arbitrary size. Rewriting Eq. (1) as

with the number of reactions *n _{r}* =

*n*+

_{g}*n*, we split the factorization into two major substeps:

_{f}- Determine the factorization array $Unc\xd7nc$ such that(19)$US=[Dng\xd7ngRng\xd7nf0(nc\u2212ng)\xd7ngE(nc\u2212ng)\xd7nf],$
where

**D**is a diagonal array with no zero entries on the diagonal,**R**can be arbitrary because it will be multiplied by*ϵ*→ 0 during the solution of the singular perturbation problem, and**E**is in upper-echelon form.Computationally,

**U**is found using integer arithmetic, partial (row) pivoting when a zero pivot element is encountered, and division of each row by its greatest common factor after each operation.If an array

**U**of form Eq. (19) can be found (i.e.,**D**is full rank), the outer solution of the singular perturbation problem can be computed by setting*ϵ*= 0 in the factored version of Eq. (18) resulting in*n*algebraic equations, in addition to_{g}*n*−_{c}*n*− rank(_{g}**E**) linear equations produced by conserved quantities.Nonfull-rank instances of

**D**indicate structural problems in the equilibrium relationships of the deposition model, such as when the modeler attempts to define independent equilibrium relationships*g*when the reaction stoichiometry requires some degree of dependence, such as in the isomerization reaction of Wei and Prater._{i}^{15}

- Rewriting Eq. (19) as$US=[Ugng\xd7ncUf(nc\u2212ng)\xd7nc][Sgnc\xd7ngSfnc\xd7nf],$the dynamic and conserved modes (
**x**and**w**, respectively) can be written as$Ufdcdt=[x\u0307w\u0307]=UfSff.$We now define the diagonal array $Jnc\xd7nc$ where

*J*_{i,}_{i}= 1 if species*c*is found in significant concentration relative to the other major species, and_{i}*J*,_{i}_{i}= 0 if species*c*is insignificant (minor)._{i}Any element of the column vector $UfJc\u0307$ that vanishes indicates a new pseudoequilibrium relationship between a subset of the finite-rate processes

*f*. Note that vanishing refers to a true zero value because of the integer arithmetic used in the factorization procedure._{i}The complete set of pseudoequilibrium relationships is computed from the row span of array $UfJ$ using a procedure akin to the Gram–Schmidt orthogonalization. This procedure is carried out upwards from the bottom row without pivoting and again, all computations are carried out using integer arithmetic. Because this process will factor out all redundant rows of $UfJ$, the maximum number of new pseudoequilibria of finite-rate processes will be extracted resulting in the minimal dynamic dimension. As with the first step, all of these operations can be represented by a single matrix, in this case $V(nc\u2212ng)\xd7(nc\u2212ng)$.

## III. RESULTS: MODELING ALUMINA ALD

To demonstrate the two-stage factorization algorithm, we consider the alumina ALD model of Remmers^{12} based on the trimethylaluminum/water precursor system. The 19 reactions for the 23 surface and gas-phase [denoted by (g)] species are listed in Table I where the finite-rate processes are shown as shaded and the true equilibrium relationships are shown with the lighter background. The species' notation is used to denote O and Al atoms with coordination numbers lower than the bulk film Al and O atoms, as discussed in detail in previous publications of the author.^{12,14} Species with primed atoms generally are considered the more reactive surface and gas-phase species. Reaction activation energies and critical complex configurations were taken from three quantum chemistry studies.^{7,8,10}

f_{0} | $Me3Al\u2032(g)+HO\u2032+3S$ | $\u21cc$ | Me_{3}AlHO |

g_{1} | Me_{3}AlHO | $\u21cc$ | Me_{3}AlHO^{‡} |

f_{2} | Me_{3}AlHO^{‡} | → | $Me2Al\u2032+Al\u2032+Al+S+MeH(g)$ |

g_{3} | $Me2Al\u2032+HO\u2032$ | $\u21cc$ | Me_{2}AlHO |

g_{4} | Me_{2}AlHO | $\u21cc$ | Me_{2}AlHO^{‡} |

f_{5} | Me_{2}AlHO^{‡} | → | $MeAl\u2032+Al\u2032+S+MeH(g)$ |

g_{6} | $MeAl\u2032+HO\u2032$ | $\u21cc$ | MeAlHO |

g_{7} | MeAlHO | $\u21cc$ | MeAlHO^{‡} |

f_{8} | MeAlHO^{‡} | → | $Al\u2032+S+MeH(g)$ |

f_{9} | $H2O\u2032(g)+2Al\u2032$ | $\u21cc$ | H_{2}O |

g_{10} | H_{2}O | $\u21cc$ | $H2O2\u2021$ |

f_{11} | $H2O2\u2021$ | → | $2HO\u2032$ + O |

f_{12} | $Me3Al\u2032(g)+Al\u2032+3S$ | $\u21cc$ | Me_{3}AlO + Al |

g_{13} | Me_{3}AlO | $\u21cc$ | $Me2Al\u2032+MeAl$ |

g_{14} | $Me2Al\u2032+Al\u2032$ | $\u21cc$ | Me_{2}Al_{2}O^{‡} |

f_{15} | Me_{2}Al_{2}O^{‡} | → | 2MeAl |

f_{16} | $H2O\u2032(g)+2MeAl$ | $\u21cc$ | MeAlH_{2}O + MeAl |

g_{17} | MeAlH_{2}O + MeAl | $\u21cc$ | Me_{2}Al_{2}H_{2}O^{‡} |

f_{18} | Me_{2}Al_{2}H_{2}O^{‡} | → | $MeAlHO+Al\u2032+S+O+MeH(g)$ |

f_{0} | $Me3Al\u2032(g)+HO\u2032+3S$ | $\u21cc$ | Me_{3}AlHO |

g_{1} | Me_{3}AlHO | $\u21cc$ | Me_{3}AlHO^{‡} |

f_{2} | Me_{3}AlHO^{‡} | → | $Me2Al\u2032+Al\u2032+Al+S+MeH(g)$ |

g_{3} | $Me2Al\u2032+HO\u2032$ | $\u21cc$ | Me_{2}AlHO |

g_{4} | Me_{2}AlHO | $\u21cc$ | Me_{2}AlHO^{‡} |

f_{5} | Me_{2}AlHO^{‡} | → | $MeAl\u2032+Al\u2032+S+MeH(g)$ |

g_{6} | $MeAl\u2032+HO\u2032$ | $\u21cc$ | MeAlHO |

g_{7} | MeAlHO | $\u21cc$ | MeAlHO^{‡} |

f_{8} | MeAlHO^{‡} | → | $Al\u2032+S+MeH(g)$ |

f_{9} | $H2O\u2032(g)+2Al\u2032$ | $\u21cc$ | H_{2}O |

g_{10} | H_{2}O | $\u21cc$ | $H2O2\u2021$ |

f_{11} | $H2O2\u2021$ | → | $2HO\u2032$ + O |

f_{12} | $Me3Al\u2032(g)+Al\u2032+3S$ | $\u21cc$ | Me_{3}AlO + Al |

g_{13} | Me_{3}AlO | $\u21cc$ | $Me2Al\u2032+MeAl$ |

g_{14} | $Me2Al\u2032+Al\u2032$ | $\u21cc$ | Me_{2}Al_{2}O^{‡} |

f_{15} | Me_{2}Al_{2}O^{‡} | → | 2MeAl |

f_{16} | $H2O\u2032(g)+2MeAl$ | $\u21cc$ | MeAlH_{2}O + MeAl |

g_{17} | MeAlH_{2}O + MeAl | $\u21cc$ | Me_{2}Al_{2}H_{2}O^{‡} |

f_{18} | Me_{2}Al_{2}H_{2}O^{‡} | → | $MeAlHO+Al\u2032+S+O+MeH(g)$ |

The reactions of Table I describe the adsorption of TMA (trimethylaluminum) ($Me3Al\u2032$) onto surface hydroxyl groups ($HO\u2032$) with net forward rate *f*_{0}, followed by the ligand exchange reactions that release methane (MeH) and ultimately produce bare aluminum surface sites $Al\u2032$ at rate *f*_{8}. Likewise, TMA can adsorb onto the surface $Al\u2032$ sites at rate *f*_{12} and then proceed to dissociate producing three MeAl surface species through *f*_{15}. Water ($H2O\u2032$) precursor molecules can adsorb on $Al\u2032$ to repopulate the surface with hydroxyl groups ($HO\u2032$), or onto MeAl sites to ultimately release MeH and leave surface $Al\u2032$. These reactions are described in detail in Remmers^{12} and the references contained therein.

Writing a dynamic material balance for each species gives 23 ODEs in time in the form of Eq. (18). Applying step one of our factorization procedure gives nine AEs in the form of *g _{i}* = 0, six conserved quantities, and eight ODEs in time that can be integrated from an initial condition of a partially hydroxylated surface initial condition to the limit-cycle solution shown in Fig. 4; the specified initial condition is first projected onto $Q$ where the solution subsequently evolves in time. This simulation corresponds to a deposition cycle taking place at

*T*= 450 K with a process recipe consisting of a 1 s TMA exposure at 2 Pa, a 5 s purge period, a 1 s water exposure at 2 Pa, followed by another 5 s purge. Therefore, each of the precursor exposure periods corresponds to 15 000 L. We note that while there are ten finite reaction rates

*f*, our factorization procedure reveals that these dynamic modes cannot be completely decoupled for this model. This simply indicates that each of the relatively slow dynamic modes will potentially have more than one time scale, dictated by the individual reaction rate values.

_{i}### A. Further dynamic dimension reduction

We compute the mean absolute value of each species (because water and TMA are consumed, $[Me3Al\u2032]$ and $[H2O\u2032]$ are nonpositive) over the first five deposition cycles and plot the results in Fig. 5. As it is to be expected, the mean concentrations of the critical complexes are all vanishingly small. Likewise, a number of energetically unfavorable surface intermediates (e.g., $MeAl\u2032$) also are found in minute concentrations. We observe the first major gap in mean concentration is the two-order of magnitude break located approximately between 10^{−1} and 10^{−3 }nm^{−2}. The major species in concentrations above this gap consist of the set of ten surface, bulk, and gas-phase species *c*_{sig}

We remind the reader that the bulk-film and gas-phase species such as $H2O\u2032$ have units nm^{−2} because these quantities denote the number of molecules consumed/produced per unit of growth surface area.

Using *c*_{sig} and carrying out steps 2(a) and 2(b) of our reaction factorization procedure, starting from the last conserved quantity and working our way up through the conserved and ODE modes to compute the span of the slow/conserved species space, we find that four new pseudoequilibrium relationships result from taking into account the split between the major and minor species:

We immediately interpret three of the new pseudoequilibrium relationships as balances between the adsorption of a gas-phase precursor and the first irreversible reaction permanently fixing it to the growth surface for the cases of TMA (20) and water [(22) and (23)]. The interpretation of Eq. (21) is more subtle: the difference *f*_{5} − *f*_{8} represents the net production rate of surface MeAlHO in the TMA adsorption/reaction sequence corresponding to rates *f*_{0} to *f*_{8}. Because MeAlHO also is produced in the final reaction at a rate *f*_{18}, *our reaction factorization has uncovered the connection between the two ALD reaction cycles* in which one cycle starts with TMA adsorbing onto a surface hydroxyl $HO\u2032$ and the second on a bare aluminum site $Al\u2032$.

Of the eight original ODE modes, four remain

Inspecting the reactions and rates *f*_{0} to *f*_{18} makes the origin of Eq. (24) immediately clear: the rate at which four-coordinate O is added to the bulk film is simply the sum of the terminal reaction rates corresponding to the two reaction sequences involving the water precursor adsorption and reaction. Subtracting twice Eq. (26) from Eq. (25) yields

where the second equality results from the equilibrium relationship (21) identified in our second factorization step. The relationship (28) also makes perfect sense, given that it simply sums the MeH producing rates. Continuing this process, we subtract Eq. (25) from Eq. (27) to find the consumption of TMA governed by

Interestingly, this results in

where *f*_{2} − *f*_{5} represents the net accumulation rate of Me_{2}AlHO produced in the reaction sequence corresponding to TMA adsorbing and reacting on a (partially) hydroxylated surface. The net rate *f*_{12} − *f*_{15} corresponds to the production of $Me2Al\u2032$ by the second reaction loop where TMA adsorbs onto bare aluminum sites. Again, the factorization procedure gives new insight into the interconnections between separate deposition reaction pathways, and can result in modeling ODEs that cannot be readily obtained from the original material balance equations.

### B. Physical interpretation of the conserved modes

Six conserved quantities are identified in the first step of the model reduction procedure; they do not change during the second step. Written in terms of the major species

The first (31) simply reflects the conservation of the Me (CH_{3}) groups; each of the reactions leaves the Me units unchanged and so the total number of Me groups must be conserved with total value *w*_{0}. The second mode (32) indicates that all of the Al incorporated into the bulk film originates in the TMA precursor. Note that this is not an elemental balance on Al and $Al\u2032$ species because of the manner in which Al generation was linked to the adsorption process in the original modeling equations. Combining Eqs. (31) and (33), we find the balance of surface sites S open to Me occupation depends on the surface concentrations of the two major Me-containing species

which, like the corresponding conserved mode for Al, simply reflects the manner in which we model the incorporation of O into the bulk film. Adding Eqs. (34)–(36) gives

which assures the conservation of total reactive sites on the deposition surface. Finally, solving Eq. (38) for [O] and substituting the result into Eq. (34), we find

which represents the conservation of H throughout all of the H-transfer reactions.

### C. Potential for lower-dimensional dynamics

An additional benefit of analyzing the conserved quantities of Sec. III B is to guide plotting the ALD system solution dynamics in one, two, or three dimensions. For example, by Eq. (38), we see that there would be no benefit to plotting the surface reaction dynamics in the ($[O],\u2009[H2O\u2032]$) plane: the dynamics in this case are restricted to the line defined by Eq. (38). Likewise, three-dimensional plots with coordinates ([MeAl], [Me_{2}AlHO], and [S]) would result in all trajectories being limited to the plane defined by Eq. (37).

Our choice is to plot the surface reaction dynamics in ($[Al\u2032],\u2009[HO\u2032],\u2009[MeAl]$)-space; we can analyze the set (31)–(36) to show no linear combination of the conserved modes defines a plane in this species subspace. Representative dynamics illustrating five ALD cycles is shown in Fig. 6. The point marked by a circle in the $([Al\u2032],[HO\u2032])$-plane indicates the initial condition $[Al\u2032](t=0)=[HO\u2032](t=0)=5.75\u2009nm\u22122$ and [MeAl](*t* = 0) = 0, a surface considered to be 50% hydroxylated. During the first TMA exposure period, both the surface Me and $Al\u2032$ site concentration increase while $[HO\u2032]$ drops. The barely visible green segment that follows shows the minor changes that occur during the purge period in this formulation of the model (reduced rates of reactions such as *f*_{8} will make the purge-period reaction dynamics more significant). As expected and as seen in Fig. 4, $[HO\u2032]$ increases dramatically during the water exposure while surface concentrations of MeAl and $Al\u2032$ plunge. Again, the purge-period changes are relatively insignificant. As we continue the ALD cycles, the reaction system settles to its limit cycle solution; what is significant in this plot is that the limit-cycle solution appears to be approximately confined to a two-dimensional, nonlinear (nonplanar) surface in this three-dimensional species space. Whether this indicates further reduction in dynamic dimension is possible will be addressed in follow-up work.

### D. Measurable species

As a final indication of the utility of our modeling and model reduction approach, consider Fig. 7 where the dynamics of the bulk film species O and Al, the cumulative consumption totals of gas-phase precursors $Me3Al\u2032$ and $H2O\u2032$, and the cumulative by-product methane MeH production totals are plotted for three full ALD cycles. We recall that precursor curves take on negative values because these values denote consumption totals of precursor molecules per nm^{2} of growth surface. After most of the surface initial state characteristics are eliminated during the first ALD cycle, we observe the regular increases in Al and O cumulative totals during the TMA and water exposure periods. These values allow us to compute the growth per cycle *gpc* as a function of cycle number (for the representative simulation we consider in this work we find *gpc* = 1.23 Å/cycle) given the density of amorphous Al_{2}O_{3}. Furthermore, we observe that O accumulates at a rate 1.5 times that of Al, indicating the correct stoichiometry of the deposited film. A similar approach that would also include the mass of adsorbed water, and chemisorbed Me and OH groups would allow the instantaneous comparison of modeling and real-time quartz-crystal microbalance or ellipsometry measurements for model validation.

We also observe the growing, negative values of $[Me3Al\u2032]$ and $[H2O\u2032]$ and the clear correspondence between these values and the dynamics of [Al] and [O], indicating a material balance for the deposited species. Finally, we observe the increasing [MeH] cumulative values; the overall molar rate of MeH production is, on average, three times that of $MeAl\u2032$ consumption, which again indicates a correct material balance. It is interesting to note that the majority of MeH production occurs during the TMA exposure, a prediction that can be compared to residual gas analysis experimental measurements.

## IV. CONCLUSION

The numerical techniques developed in this work were motivated by the current lack of a coherent approach to integrating reaction path segments from a range of computational chemistry studies into a viable CVD or ALD reaction network and to distill that RN to its lowest dynamical dimension. The ultimate objective of this line of work is to create a constructive set of numerical tools that produce physically correct deposition reaction networks and the corresponding predictive computational models intended to be useful for high-throughput thin-film process design and optimization.

Our approach consists of two major steps: (1) to identify through a reaction factorization procedure redundant dynamic modes and to simultaneously formulate a proper singular perturbation problem for systems with reaction equilibrium constraints, and (2) define major/minor species to identify finite-rate processes at pseudoequilibrium. An important consequence of the first step is that if the equilibrium reactions cannot be decoupled, the problem is identified as nonphysical and so must be reformulated. Additionally, physically meaningful conserved quantities are identified in the first step of the procedure, as well as in the second step where the new algebraic relationships can be attributed to pseudoequilibria between adsorption and reaction processes as well as interconnections between ALD reaction network loops that would not be easily discerned without this approach.

The numerical methods developed in this work were applied to an alumina ALD process based on the TMA/water precursor system. The RN for this process consisted of 19 reactions and 23 surface and gas-phase species, a model that ultimately was reduced to four dynamical dimensions. For the representative operating conditions of *T* = 450 K and 15 000 L exposures of each precursor both the reduced and original models predicted saturating ALD growth dynamics, a *gpc* = 1.23 Å/cycle, and the correct film alumina film composition, results consistent with experimental observations. While the original and reduced models produced essentially identical simulation results, validating our reduction approach, a comparison to experimentally observed reaction dynamics was beyond the scope of this work. However, an advantage of the analysis and reduced model development procedure described is its ability to predict real-time surface species composition, film mass dynamics, and reactant by-product evolution rates, all potentially valuable for model validation.

Because of the goal of producing models useful for ALD process design and optimization, current computational efforts must focus on automated generation of the reduced models in a suitably self-contained form. Additional automation of the model analysis, such as interpretation of the conserved modes and pseudoequilibrium relationships, also would be of real benefit to the model builder. Finally, extension to examples beyond alumina ALD is underway, particularly to ALD processes producing crystalline films.

## ACKNOWLEDGMENTS

The author gratefully acknowledges the support of the U.S. National Science Foundation through Grant Nos. CBET1160132 and CBET1438375.