The InsTITuTe for sysTems research Isr develops, applies and teaches advanced methodologies of design and analysis to solve complex, hierarchical, heterogeneous and dynamic prob- lems of engineering technology and systems for industry and government. Isr is a permanent institute of the university of maryland, within the a. James clark school of engineering. It is a graduated national science foundation engineering research center. www.isr.umd.edu Reaction Path Analysis for Atomic Layer Deposition Systems Raymond Adomaitis Isr TechnIcal rePorT 2018-02 REACTION PATH ANALYSIS FOR ATOMIC LAYER DEPOSITION PROCESSES Raymond A. Adomaitis⇤1 1Department of Chemical & Biomolecular Engineering, Institute for Systems Research, University of Maryland, College Park, MD 20742, USA Abstract In this paper, we examine the mathematical structure of thin-film deposition process reaction kinetics models with the goal of determining whether a reaction network can guarantee the self-limiting and stable growth inherent in true atomic layer deposition systems. This analysis is based on identifying reaction invariants and interpreting the chemical significance of these conserved modes. A species-reaction graph approach is introduced to aid in distinguishing “proper” from problematic ALD reaction networks. Keywords Atomic layer deposition, reaction invariants, singularly perturbed system, species-reaction graph. Introduction In atomic layer deposition (ALD) processes, the growth surface is exposed to cycles of alternating gas- phase precursors to produce thin solid films with atomic- level thickness control. Because of the ability of ALD to deposit an increasingly wide range of elements and compounds over topographically varying surfaces with film morphologies ranging from amorphous to crystalline (Miikkulainen et al., 2013), ALD is emerging as a critical manufacturing technology for energy storage and con- version, nanoelectronics, and biomedical applications (George, 2010). Unlike its chemical vapor deposition (CVD) coun- terpart, ALD processes have no steady mode of opera- tion, and so process optimization requires kinetics mod- els of the deposition reaction network (RN). Significant progress has been made in modeling ALD surface pro- cess from a first-principles perspective (Elliott, 2012), but studies of complete ALD RN are rare because 1. Many ALD kinetics studies are limited to a por- tion of the full RN (see, e.g., Travis and Adomaitis (2013)) resulting in fragmented reaction mecha- nisms studies; 2. There are competing reaction paths to a prod- uct species involving widely ranging, multiple time ⇤To whom all correspondence should be addressed ado- maiti@umd.edu scales (Delabie et al., 2012); 3. The mechanistic origins of self-limiting and steady cyclic-growth processes remains open to debate (Puurunen, 2005). The objective of this paper is to develop the analysis tools necessary to assess whether we have an “proper” ALD reaction network before investing in the substantial e↵ort of determining reaction rates. Elements of an ALD reaction process Let us consider an archetype ALD process consist- ing of the following overall reaction between a metal- containing precursor ML2 and water to produce a metal- oxide film MO and gas-phase by-product HL: ML2(g) + H2O(g) ! MO(b) + 2HL(g) This reaction can represent, for example, the ALD of ZnO from diethyl zinc and water precursors (Gao et al., 2016). The elementary reactions and net-forward re- action rates for each of the half-reactions are given in Table 1. The first three reactions of Table 1 correspond to the metal precursor half-reaction. This simple sequence of reactions begins with the reversible adsorption of ML2 onto the O of a surface hydroxyl group HO with net forward rate f0 to produce the surface adduct HML2. Because this adsorption reaction renders the O inert to Reaction Net rate s1 m2 ML2(g)+2S+HO ⌦ HML2+O(b) f0 HML2 ⌦ HML‡2 (1/✏)g1 HML‡2 ! HL(g)+S+ML f1 H2O(g)+ML ⌦ H2OL+M(b) f2 H2OL ⌦ H2OL‡ (1/✏)g2 H2OL‡ ! HL(g)+S+HO f3 Table 1. Archetype ALD process reactions and rates. all subsequent reactions (except, of course, the desorp- tion of ML2), we denote the incorporation of O into the bulk film by the production of species O(b). Addition- ally two surface sites S are consumed by this reaction; this fictitious species accounts for the area of reaction surface that is sterically hindered by the two metal pre- cursor ligands L. We note that gas-phase and bulk-film species are explicitly noted as such by (g) and (b), re- spectively, while all others are surface species. The adsorbed adduct HML2 can undergo a (1-2) H- transfer reaction (Delabie et al., 2012) by forming the critical complex HML‡2 involving the H of the hydroxyl group onto which the precursor adsorbed. We note that while conventional transition-state theory (CTST) dic- tates this to be an equilibrium process (Laidler, 1987), we write the net-forward reaction rate as the finite- rate process (1/✏)g0 using relaxation time constant ✏ for the purpose of correctly formulating the species bal- ances that follow. The transition state HML‡2 then can eliminate by-product HL(g) and liberate an adsorption site S through finite-rate process f1, which leaves the permanently-bonded surface species ML. For this study, we consider this reaction irreversible because of sucient reactor exhaust rate to e↵ectively remove all HL. The reactions corresponding to the water exposure mirror those of the metal-precursor half-cycle: water adsorbs onto the reactive metal surface species ML to form adduct H2OL, resulting in the incorporation of metal M(b) into the bulk film. However, because there is no change in surface ligand concentration [L], there is no corresponding consumption of S. Again, the surface adduct can undergo a H-transfer reaction by first form- ing critical complex H2OL‡ in an equilibrium reaction governed by g1, and then can undergo an irreversible reaction ejecting another HL(g), freeing one surface site S and leaving a surface HO group by finite-rate process f3. Before proceeding, we note that reaction rates of Ta- ble 1 can be generated from experimental measurements or quantum chemical computations coupled to conven- tional transition-state theory (Laidler, 1987); the anal- ysis that follows is entirely independent of values of re- action rates fi and the nature of the gi. Atoms, species, phases, and balances The twelve species (including surface sites S) of Table 1 can be collected into chemical species set S: S =HML2,HML‡2,H2OL,H2OL‡,ML2(g), S,HO,O(b),HL(g),ML,H2O(g),M(b) . (1) Clearly, three phases exist for this reaction system: the gas, surface, and bulk film phases P = {0 (gas), 1 (growth surface), 2 (film)} (2) where 0 corresponds to the gas volume in nm3 and 1 the reaction surface area in nm2. While these quanti- ties are necessary to define the molar species balances, we are free to choose 2 to either represent the total film volume or the film surface area - the latter case is useful when we wish to represent the number of bulk species (O and M) incorporated per unit area of the growth sur- face. From the reactions listed in Table 1 we also can extract a set of four “elements” E = {M,O,L,H} (3) where the ligand L is included in E because it remains untransformed by any of the proposed surface reactions. As such, we note that the notation used for the chemical species in S can be thought of as first step towards mov- ing to a highly simplified form of the SMILES notation (Weininger, 1988). The thermodynamic system we study is a di↵erential volume of constant size 0 that is perfectly uniform in each of its phases and is closed to the environment; at this time we place no restrictions on the reaction sur- face other than 1 > 0 at initial conditions t = 0. With the species S (1), phases P (2), and the reactions, stoi- chiometry, and reaction rates of Table 1 in hand, we can write the twelve species di↵erential equation balance as dm dt = 1 ✏ P " g0 g1 # + 1Q 26664 f0 f1 f2 f3 37775 (4) with m being the molar amounts (not concentrations) of each species in ns ⇥ 1 array arranged in the ordering of (1), subject to specified initial conditions (Remmers et al., 2015): n(t = 0) = nAo , n(t = ⌧ A + ⌧AP ) = nBo (5) at the start of the ML2 (at t = 0) and water doses (t = ⌧A + ⌧AP ) respectively, where ⌧A is the length of the ML2 exposure and ⌧AP the post-ML2 purge period. The stoichiometric arrays Pns⇥ng and Qns⇥nf are: P = 2666666666666666666664 1 0 1 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3777777777777777777775 , Q = 2666666666666666666664 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 1 0 0 0 2 1 0 1 1 0 0 1 1 0 0 0 0 1 0 1 0 1 1 0 0 0 1 0 0 0 1 0 3777777777777777777775 (6) Given the underlying assumption of CTST that the re- actions producing critical complexes HML‡ and H2OL‡ are in equilibrium, the true solution to (4) is found my multiplying the di↵erential equations through by ✏ and then taking the limit ✏ ! 0. Because the first four bal- ances have nonzero entries in P, this operation results in the loss of two equations, making it impossible to solve (4). Despite the simplicity of our archetype ALD model, computing its solution is much more complicated than one might suspect because (4) constitutes a singularly perturbed system (STS) in non-standard form (Daou- tidis, 2015). A “proper” ALD RN This mechanistically simple but mathematically non- trivial model opens the question of how one defines a “proper” ALD process model. Focusing only on the in- trinsic deposition kinetics, we pose the following as a subset of questions that must be asked of an RN model structure before significant e↵ort is expended on identi- fying reaction rates and before precursor and by-product gas-phase transport phenomena modeling elements are incorporated to complete the deposition system descrip- tion: 1. Can the reaction process time scales be separated even at the coarsest level, e.g., between equilibrium and finite-rate processes when ✏! 0? 2. Will the model reduction process used to analyze the RN model indicate whether it will be possible to measure each finite-rate reaction process indepen- dently, or will additional information beyond the time-rate of change of species in S be required to determine the reaction rate values? 3. Will the deposited film have the correct stoichiom- etry regardless of the reaction-rate values? 4. Is the overall RN balanced, i.e., are all elemental balances satisfied for all time including the system’s original state, as well as during the transitions be- tween precursor doses and the purge periods even in the limit of infinitely fast transitions? 5. Is the deposition surface stable, e.g., will the reac- tion surface area have a positive and bounded value for any number of ALD deposition cycles? 6. Are self-saturating conditions guaranteed to ex- ist, and can the mechanism be unambiguously ex- tracted from the RN without information regarding the reaction rates? 7. Will the modes identified as being dynamically re- dundant have a physical meaning in the context of ALD, and can this meaning be identified as part of the reduction process? Invariant analysis of the archetype ALD process Because our system is closed and a balance for every species in the ALD system is provided, we should ex- pect (4) to contain redundant dynamic modes because elements – and potentially other reaction quantities – must be conserved. Fortunately, these modes can be identified concurrently with the transformation of (4) to a STS in standard form through a reaction factorization (diagonalization) process (Adomaitis, 2016). Defining the arrays R = [P,Q]ns⇥nr , h = " g f #ns⇥1 with nr = ng +nf , we decouple the reactions through a Gauss-Jordan elimination procedure (Adomaitis, 2016; Remmers et al., 2015; Rodrigues et al., 2015) to as nearly as possible produce TR = " Inr⇥nr 0(nsnr)⇥nr # (7) where array T is the matrix equivalent of the diagonal- ization procedure. The objective of “nearly as possible” actually is im- portant from a chemical kinetics point of view: under some circumstances (such as with our archetype ALD process), it is possible to achieve the transformation (7) exactly. This indicates the independence of the re- actions in this reaction network (RN). However, many RN can feature elementary reaction sequences that form competing paths to the same chemical species product - these situations constitute, of course, perfectly legiti- mate chemical RN, and result in the inability to satisfy (7) exactly. Under these circumstances, it is sucient to reduce (4-6) to upper-echelon form which transforms the STS to standard form, eliminates redundant dynamic modes, but does not decouple the finite-rate reaction terms (Remmers et al., 2015). The reaction factorization procedure can be carried out to completion for the archetype ALD RN (4-6) us- ing integer arithmetic to avoid any numerical ambiguity or loss of numerical precision resulting in a diagonal- ized system equivalent to (7). Application of the di- agonalization procedure produces ng = 2 independent algebraic relationships corresponding to the equilibrium reactions, nf = 4 ordinary di↵erential equations in time corresponding to four dynamically decoupled states xi, i = 0, 1, 2, 3, and nc = ns nr = 6 conserved quantities wi, i = 0, . . . , 5. The transformed system is given in Table 2. The di↵erential-algebraic equation (DAE) system of Table 2 now is in the form that can numerically solved using standard DAE solvers (e.g., an implicit- Euler scheme works very well for this relatively low- dimensional system for reaction rates fi that do not span a wide range of timescales) provided the initial conditions marking the onset of the metal- and water- precursor doses (5) can be projected onto g0 = 0 and g1 = 0 (these functions are linear for the archetype ALD process, and so satisfying this condition is trivial); de- tails regarding these numerical issues can be found in Remmers et al. (2015). Given the decoupled ALD reaction model, we return to the list of criteria posed in Section for a “proper” ALD process model to examine which questions have been resolved and which remain open: • The reaction diagonalization (factorization) pro- cedure unambiguously determines if the pseudo- equilibrium relationships gi = 0 can be solved in- g0 = 0 g1 = 0 d dt (ML2) = dx0 dt = f0 d dt (HML2 HML‡2 ML2) = dx1 dt = f1 d dt (HML2 +HML ‡ 2 +H2OL +H2OL ‡ ML2 + S) = dx2 dt = f2 d dt (HML2 +HML ‡ 2 ML2 + S) = dx3 dt = f3 HML2 HML‡2 S +HO = w0 ML2 +O = w1 2ML2 S +HL = w2 2HML2 + 2HML ‡ 2 +H2OL +H2OL ‡ + S +ML = w3 HML2 +HML ‡ 2 +H2OL +H2OL ‡ ML2 + S +H2O = w4 HML2 HML‡2 H2OL H2OL‡ +ML2 S +M = w5 Table 2. The archetype ALD reaction network model transformed to a singularly perturbed system in standard form, resulting in two algebraic equations, four ordinary di↵erential equations in time, and six conserved modes. dependently at all times during the simulation, in- cluding the initial conditions, confirming condition (1) of the list. • The new states xi, i = 0, . . . , 3 are linear combi- nations of the species molar quantities; these are known as reaction variants (Asbjørnsen, 1972; Ro- drigues et al., 2015; Zhao et al., 2016). The reaction variants determine the deposition system’s minimal dynamic dimension and whether those rates can be experimentally and independently measured; for the archetype ALD process, the rates can be mea- sured independently and so condition (2) of the list is satisfied. • The six modes wi, i = 0, . . . , 5 are known as the reaction invariants. These represent combinations of molar species quantities that remain invariant in time. Examining the wi defined in Table 2, how- ever, reveals little in terms of the physical meaning of these quantities. Therefore, additional analysis is required to address issues (3-7). Species-Reaction graphs for extracting and in- terpreting invariants The Species-Reaction (SR) graph was introduced by Feinberg and coworkers (Craciun and Feinberg, 2006) to facilitate analysis of chemical reaction networks to deter- mine the potential for multiple equilibria under isother- mal conditions. An SR graph is constructed from a se- quence of chemical reactions according to the following format: 1. Chemical species and reactions form the nodes of the graph; species are denoted with circular nodes, reactions with square nodes. In our work, we dis- tinguish equilibrium from finite-rate reactions using di↵erent node colors (blue and yellow, respectively), although the distinction has no bearing on the iden- tification and interpretation of reaction invariants. 2. Species and reactions are connected by edges de- fined by the stoichiometry of the reactions; as such, species can only be connected to reactions so reaction-reaction and species-species edges are not possible. 3. Reaction stoichiometric coecients are used to la- bel the edges. In our work, negative stoichiomet- ric coecient values denote reactants and positive products – the sign notation can be reversed with- out a↵ecting the results of our analysis. Translation of the archetype ALD process of Table 1 (it it important to stress that the graph corresponds to the RN prior to our factorization procedure) results in the graph shown in Fig. 1. With this graph, one can visually trace the path of reaction sequences through the RN from gas-phase precursor (e.g., ML2), through the surface reactions, to the ultimate destination of the metal atom M in the bulk film. Likewise, it also is pos- sible to trace the multiple reaction paths that lead to the gas-phase by-product HL. Terminal to terminal species, linear graph Given the complexity of the archetype ALD process SR graph (Fig. 1), we now turn to examining its distinct Figure 1. ALD archetype system SR graph. ML2 and H2O represent the gas-phase precursors, HL the gas-phase by-product, O and M are the bulk film com- ponents, and all other components constitute surface species. Critical complexes are denoted by a “+” sux. Equilibrium reactions are shaded blue, while finite-rate precesses are shown in yellow. subgraphs to determine their connections to the reac- tion invariants of the complete RN. To start, consider a simple dissociative adsorption process where gas-phase dimer species D(g) adsorbs dissociatively and reversibly onto a reaction surface; the physisorbed adatoms A then can undergo an irreversible incorporation into the film consisting of species B(b). The two reactions and three species are written as: |⌫R0|D(g)⌦ ⌫P0A with net forward rate f0 (8) |⌫R1|A! ⌫P1B with rate f1 (9) with ⌫R0 = ⌫R1 = 1, ⌫P1 = 1, and ⌫P0 = 2. These species, reactions, and stoichiometric coecients are shown as an SR graph in Fig. 2. Because all of the ter- minal nodes of this SR graph are chemical species and none is a rate process, the system is closed and so we expect at least one invariant quantity; this is easily ver- ified by writing the three di↵erential equation balances for species D(g), A, and B(b) and then performing the reaction diagonalization procedure. Alternatively, we can trace a path through the SR graph from the D to B nodes with molar quantities of D, A, and B to find D + |⌫R0| |⌫P0|  A+ |⌫R1| |⌫P1|B = constant or 2D +A+B = w0 (10) Physically, we interpret (10) as being equivalent to the conservation of the element deposited in the film by species B(b), where w0 corresponds to the sum of the moles of the element in each of the species, starting from the precursor D(g) and ending with B(b). νR0 νP0 νR1 νP1 Figure 2. Linear species-reaction graph corresponding to the dissociative adsorption reactions (8-9). Reaction branches Now consider the SR graph of the single reaction H2 +AB! AH+ BH (11) shown in Fig. 3. As with the previous reaction system, all terminal nodes of the SR graph corresponding to (11) are species, and so Fig. 3 represents a closed system. There are a total of six paths we can take from one terminal species to another; diagonalization of the four species balances reveals three conserved modes and only one independent species balance. Four of the paths pass through the reaction: a) H2 + AH = constant b) H2 + BH = constant c) AB + AH = constant d) AB + BH = (b) + (c) - (a) and so the fourth path (d) is clearly seen as a linear combination of the first three. Two paths bypass the reaction and so con- stitute paths through each reaction complex: H2 - AB = (a) - (c) BH - AH = (b) - (a) which are both linear combinations of the paths (a-d) passing through the reactions. Figure 3. Reaction branch SR graph for (11). At this point, we make one observation relevant to the interpretation of the reaction invariants of Ta- ble 2: that the negative quantities originate from reac- tion branches in the archetype ALD process SR graph representing paths through reaction complexes. To gen- erate physically meaningful reaction invariants, we will only generate reaction invariants corresponding to paths passing through reactions (so the stoichiometric coe- cients must change sign between the incoming and out- going edges), and to combinations of those modes result- ing from sums of invariants. For example, because (11) must produce three independent reaction invariants, we can add paths (a) and (b) and list (c) and (d) as-is to define 2H2 + AH + BH = w0 AB + AH = w1 AB + BH = w2 which correspond to atomistic balances of H, A, and B summing to w0, w1, and w2, respectively. We observe that this path-following procedure to generate invari- ants in some ways can be seen as a “logical or” for the evolutionary paths of chemical species in the RN - that certain species potentially participate in one “or” more conserved relationships, and the enumeration of those paths takes place in a manner independent of other re- action paths corresponding to conserved quantities. A final connection between well-established reaction stoichiometry relationships and reaction invariants can be observed by creating the atomic balance array A: H2 AB AH BH A 0 1 1 0 (w1) B 0 1 0 1 (w2) H 2 0 1 1 (w0) where the nullity of A corresponds to the number of columns of A - the rank of A, which in this case has the value 1. This means the null space or kernel of A can be found as the one-dimensional vector [1,1, 1, 1]T which is, of course, a vector of the stoichiometric coe- cients of (3). Species branches Next consider the RN of Fig. 4, where the primary di↵erence relative to the previous case is that this RN branches through a species node as opposed to a reaction as in Fig. 3. As an RN, this represents the (reversible) conversion of species A to B, the latter of which can be converted to either species C or D. The selectivity to these terminal species is determined by the relative rates f0 and f1. Figure 4. Species branch graph corresponding to the re- action network A ⌦ B, B ! C, B ! D. Formulating the species balances and diagonalizing, we find only one reaction invariant. Conceptually, we can split the total number of A and B which ultimately are converted to C and D by defining A = A(0) + A(1) and B = B(0) + B(1) where superscripts (0) and (1) denote that subtotals of A and B destined to become C and D, respectively. Following each RN path from A to either C or D gives A(0) +B(0) + C = constant (12) A(1) +B(1) +D = constant (13) therefore A+B + C +D = w0 (14) It is important to stress that because species molar amounts A(0), A(1), B(0) and B(1) are fictitious quanti- ties, (12) and (13) are not physically realizable invari- ants; equation (14) is the only true reaction invariant of this RN. As such, because all paths passing through a species branch must be summed, this type of branch serves as a “logical AND” for paths through a RN. Reaction cycles The final RN structure we examine is the reaction cycle. Consider the sequence of reversible decomposition reactions A ↵ 2B ↵ 4C shown in Fig. 5. We close this sequence of reactions into a cycle by allowing 4C ↵ A; this is, of course, highly contrived, but that is intentional so as to make clear what role the reaction stoichiometry plays in the reaction invariant associated with this RN. We note, however, that the structure of this cyclic reaction is exactly the same as that of the isomerization reaction studied byWei and Prater (1962). As with the species-branching case, we split the num- ber of one species into two artificial subtotals, such as A = A(0) +A(1) which breaks the cycle into a terminal- to-terminal linear reaction sequence. The procedure to account for the reaction stoichiometry used in (10) then Figure 5. Reaction cycle graph for A ↵ 2B ↵ 4C. produces the single reaction invariant A(0) + 1 4  C + 2 1  B + 2 1 A(1) = constant 4A+ 2B + C = w0 Interpretation of this result makes physical sense: if species A corresponds to the element making up species B and C, w0 represents the initial and constant number of atoms of A in the RN. Archetype ALD process: SR graph analysis Returning to the archetype ALD process described by the reactions of Table 1 and shown as a SR graph in Fig. 1, we now present an alternative approach to determining the six expected reaction invariants and will examine their physical meaning relative to those listed in Table 2. As seen in Fig. 6, there are two non- branching paths (shaded in blue) in this RN: one orig- inates from the metal gas-phase precursor ML2(g) and travels through the SR to the bulk metal M(b), while the other path originates from the O-containing precursor H2O(g) and follows a di↵erent path to bulk film O(b). These are shown as the top two SR graphs of Fig. 6; the corresponding reaction invariants are given as w0 and w1 in Table 3. Two branching paths can be seen in the lower graphs of Fig. 6. As with the M and O conservation modes, the two new paths have as one terminus either the ML2(g) or the H2O(g) precursor. However, each of these cases possesses a reaction branch point, where both branches lead to the HL(g) by-product (these paths are di↵eren- tiated by color in Fig. 6). This ultimate confluence of reaction paths can be treated either as a species branch or as a cycle; in either case, the reaction invariant is Figure 6. Archetype ALD SR graph illustrating reaction paths corresponding to M (top left), O (top right), H (bottom left), and L (bottom right) conservation. resolved by the definition of the HL(g) molar subsets HL = HL(0) +HL(1), allowing for the identification of two invariants for each case which are then added to eliminate the HL(0) and HL(1). This analysis results in the L and H conservation modes, defining the third and forth invariants (w2 and w3) of Table 3. The surface SR subgraph Having identified four physically meaningful and lin- early independent reaction invariants, two remain to be found. To aid in this identification process and clarify the physical interpretation of these modes, we remove the gas- and bulk-phase species from the RN SR graph, leaving only the surface species and reactions. The sim- plified RN SR graph is shown in Figs. 7 and 8. Inspec- tion of the SR graph limited only to surfaces processes reveals three interconnected cycles that define the two remaining reaction invariants. Two of these cycles are shown in Fig. 7, where the overlap between the cycles consists of the S - f0 - HML2 - g0 - HML ‡ 2 - f1 path. Both cycles can be interpreted as closed loops of reaction processes, each consuming and then producing one unit of surface area corresponding to the size of ligand L. Thus, one cycle (depicted in red) consumes one surface site S through the formation of adsorbed HML2 which then is transformed to critical complex HML‡2 and then consumed to return one site S as part of this H-transfer reaction. Likewise, the same sequence initiates the second cy- cle (shown in blue) corresponding to the L group not involved in the first cycle. Therefore, this cycle proceeds through the ligand-exchange path involving reactions related to the production and consumption of H2OL, which ultimately returns the second S consumed during the ML2(g) adsorption. As such, this overall process can be described as a species-branch process, where splitting species S to define the sum S = S(0)+S(1) decouples the cycles and subsequently determines the reaction invari- ML2 +HML2 +HML ‡ 2 +ML+M = w0 (M conservation) H2O +H2OL+H2OL ‡ +HO +O = w1 (O conservation) H2O +H2OL+H2OL ‡ +HL(0) + H2O +H2OL+H2OL ‡ +HL(1) +HO +HML2 +HML ‡ 2 = w2 (H conservation) ML2 +HML2 +HML ‡ 2 +HL (0) + ML2 +HML2 +HML ‡ 2 +HL (1) +ML+H2OL+H2OL ‡ = w3 (L conservation) HML2 +HML ‡ 2 + S (0) + HML2 +HML ‡ 2 + S (1) +ML+H2OL+H2OL ‡ = w4 (reaction area conservation) HML2 +HML ‡ 2 +ML+H2OL+H2OL ‡ +HO = w5 (reaction site conservation) Table 3. Conserved modes for the archetype ALD reaction system. The sums HL(0) +HL(1) and S(0) + S(1) are replaced with HL and S, respectively, in the true invariants. Figure 7. Surface-phase reaction surface area conserva- tion for the archetype ALD process. ant w4 of Table 3. Defining S = S(0) + S(1) also splits the edge between species S and reaction f0 in the SR graph into two, each with a stoichiometric coecient of 1. Identification of self-limiting ALD behavior During an exposure to the metal-containing precursor ML2(g) , reaction rate f2 = 0, a condition correspond- ing to zero H2O adsorption rate. This e↵ectively breaks the second cycle described above, resulting in a terminal species ML. This corresponds to the ALD growth sur- face for a saturating ML2(g) dose and demonstrates the self-limiting nature of a true ALD process. An alterna- tive way of understanding this behavior is to compare it a chemical vapor deposition (CVD) reactor operating at Figure 8. Surface-phase reactive site conservation for the archetype ALD process. steady state. Under this condition, one or more contin- uous cycles must exist containing species S to maintain open adsorption sites for the CVD precursors. Surface site conservation A third cycle limited to surface species and reactions can be identified that is linearly independent of the two cycles described above. As seen in Fig. 8, this clearly defined cycle does not involve S but does contain the surface hydroxyl HO. By following the reactions in this cycle, it is clear that it represents the final invariant: the conservation of reaction sites and conserved quantity w5 of Table 3. The importance of this reaction invariant is that it guarantees the reaction surface remains bounded - that it does not growth indefinitely or vanish, halting the reaction process. As with the the reaction invariant signaling self-limiting ALD behavior, this cycle also is broken when f0 = 0 and/or f2 = 0 corresponding to the individual precursor doses and purge periods. Concluding remarks The primary objective of this work was to define what constitutes a “proper” atomic layer deposition (ALD) reaction kinetics model through the physical in- terpretation of the reaction invariants. By constructing species-reaction (SR) graphs of the ALD reaction net- work (RN) and by developing a set of rules for extract- ing reaction invariants from the SR graphs, six reaction invariants were identified for our archetype ALD pro- cess. Four invariants were found to be attributed to the species elemental balances and the remaining two to re- action surface area and species conservation modes; the latter were interpreted as signals of “proper” ALDmodel behavior, both in terms of self-limiting ALD growth and stability of the growth surface. This scope of this study was limited to reaction in- variants; interpretation of the reaction variants in the context of the SR graphs is underway. The ALD RN models of this study correspond to closed systems and so do not account for the dynamics associated with the transport of reactants and gas-phase reaction by- products into and from the ALD reactor vessel. While gas-phase transport will not a↵ect the invariants associ- ated with the surface processes, analysis of the open re- action system will be necessary to understand the com- plete ALD RN picture. Acknowledgments The author gratefully acknowledges the support of the US National Science Foundation through grants CBET1160132 and CBET1438375. References Adomaitis, R. A., Dynamic dimension reduction for thin- film deposition reaction network models, Proc, IFAC DYCOPS-CAB 2016, Trondheim, Norway 448-453 (2016). Asbjørnsen O. A., Reaction invariants in the control of con- tinuous chemical reactors, Chem. Engng Sci. 27 709-717 (1972). Craciun, G. and M. Feinberg, Multiple equilibria in com- plex chemical reaction networks: II. The species-reaction graph, SIAM J. Appl. Math., 66 1321-1338 (2006). Daoutidis, P., DAEs in model reduction of chemical pro- cesses: An overview, Surveys in Di↵erential-Algebraic Equations II A. Ilchmann and T. Reis (eds), Springer In- ternational Publishing (2015). Delabie, A., S. Sioncke, J. Rip, S. Van Elshocht, G. Pour- tois, M. Mueller, B. Beckho↵, and K. Pierloot, Reaction mechanisms for atomic layer deposition of aluminum ox- ide on semiconductor substrates, J. Vac. Sci. Technol. A 30 01A127-1 (2012). Elliott, S. D., Atomic-scale simulation of ALD chemistry, Semicond. Sci. Technol. 27 074008 (2012). Gao, Z., F. Wu, Y. Myung, R. Fei, R. Kanjolia, L. Yang, and P. Banerjee, Standing and sitting adlayers in atomic layer deposition of ZnO, J. Vac. Sci. Technol. A 34 01A143-1– 10 (2016). George, S. M., Atomic layer deposition: An overview, Chem. Rev. 110 111-131 (2010). Laidler, K. J., Chemical Kinetics, Third edition Harper & Row, (1987). Miikkulainen, V., M. Leskela¨, M. Ritala and R. L. Puurunen, Crystallinity of inorganic films grown by atomic layer de- position: Overview and general trends, J. Appl. Phys. 113 021301-1–101 (2013). Puurunen, R. L., Surface chemistry of atomic layer deposi- tion: A case study for the trimethylaluminum/water pro- cess, Appl. Phys. Rev. 97 121301-52 (2005). Remmers, E. M., C. D. Travis, and R. A. Adomaitis, Re- action factorization for the dynamic analysis of atomic layer deposition kinetics, Chem. Engng Sci., 127 374-391 (2015). Rodrigues, D., S. Srinivasan, J. Billeter, D. Bonvin, Variant and invariant states for chemical reaction systems, Com- puters and Chemical Engineering, 73 23-33 (2015). Travis, C., D. and R. A. Adomaitis, Modeling alu- mina atomic layer deposition reaction kinetics during the trimethylaluminum exposure, Theo. Chem. Accounts 133 1414-1–11 (2013). Wei, J. and C. D. Prater, The structure and analysis of com- plex reaction systems, Adv. Catalysis 13 203-392 (1962). Weininger, D, SMILES, a chemical language and informa- tion system. 1. Introduction to methodology and encoding rules, J. Chem. Inf. Comput. Sci. 28 31-36 (1988). Zhao, Z., J. M. Wassick, J. Ferrio, and B. E. Ydstie, Reaction variants and invariants based observer and controller de- sign for CSTRs, Proc, IFAC DYCOPS-CAB 2016, Trond- heim, Norway 1091-1096 (2016).