Dynamic Modeling for the Design and Cyclic Operation of an Atomic Layer Deposition ( ALD ) Reactor

A laboratory-scale atomic layer deposition (ALD) reactor system model is derived for alumina deposition using trimethylaluminum and water as precursors. Model components describing the precursor thermophysical properties, reactor-scale gas-phase dynamics and surface reaction kinetics derived from absolute reaction rate theory are integrated to simulate the complete reactor system. Limit-cycle solutions defining continuous cyclic ALD reactor operation are computed with a fixed point algorithm based on collocation discretization in time, resulting in an unambiguous definition of film growth-per-cycle (gpc). A key finding of this study is that unintended chemical vapor deposition conditions can mask regions of operation that would otherwise correspond to ideal saturating ALD operation. The use of the simulator for assisting in process design decisions is presented.


Introduction
Atomic layer deposition (ALD) is a thin-film manufacturing process in which the growth surface is exposed to an alternating sequence of gas-phase chemical precursor species separated by purge periods to prevent gas-phase reactions [1,2].ALD is characterized by self-limiting heterogeneous reactions between the gas-phase precursor species and surface-bound species, which, when allowed sufficient conditions to reach saturation, results in highly conformal thin films, on both planar and non-planar geometries, with atomic level control of film deposition [1,3].With advances in current technologies alongside a growing body of knowledge on the ALD process, ALD functions have expanded to accommodate a wide range of applications, including photovoltaics [4], energy devices [5][6][7], nanofabrication [8], environmental issues [6] and even medical devices and biological systems [9].
A simplified view of the ALD binary reaction sequence for growing metal-oxide compound, "MO", from the metal precursor (denoted M with a red dot) and the oxygen precursor (denoted O with a blue dot) is schematically illustrated in Figure 1.Generally, during the first exposure period (Exposure A), the metal precursor adsorbs onto an oxygenated substrate, undergoing a ligand-exchange reaction with the surface-bound oxygen species and, thereby, becoming permanently bound to the growth surface.After a sufficient purge period leaving no reactive species in the gas phase, the oxygen precursor (e.g., H 2 O, O 2 , O 3 ) is introduced during Exposure B, initiating a subsequent reaction, which proceeds analogously to the reaction in Exposure A, whereby the oxygen precursor adsorbs onto the surface, undergoes a ligand-exchange reaction with the surface-bound metal species and, consequently, becomes permanently bound to the growth surface [3].play a crucial role in the deposition process [2,10].Careful chemical and mathematical analysis is required to resolve the multiple time scales present in this process to identify the rate-limiting steps.
The focus of this paper will be on the development of a dynamic ALD process model (Figure 2) based on a laboratory-scale reactor system currently under construction.The primary contributions of this paper are its use of physically-based reaction kinetics models derived from transition state theory [11] and limit-cycle calculations describing the steady operation of these reactor systems.

. Review of ALD Models
A range of ALD models has been developed at various levels of detail and theory [12] to accommodate specific modeling objectives.Of these, atomic-scale models have produced insightful results by employing ab initio methods, such as density functional theory (DFT), to investigate reaction mechanisms, energetically favored reaction pathways, proton transfer mechanisms, structures of surface reaction intermediate species and reaction energetics corresponding to reactions along each stage of ALD growth (see [13] for a comprehensive review of atomic-scale ALD models).Results of such studies have been combined with kinetic Monte Carlo (kMC) simulations to quantify nucleation [14] and growth [15,16] kinetics.Empirical growth models have also been developed to describe ALD kinetics during the initial stages of nucleation [17,18] and island growth [19][20][21].Additionally, surface reaction models have been developed to predict growth kinetics and the effects of operating parameters (e.g., temperature, pressure, precursor exposure time) during undersaturating conditions and leading to self-limiting ALD growth [11,[22][23][24][25].
A number of detailed ALD reactor simulations have been reported in the literature, generally comprised of transport equations coupled with a surface reaction model to describe the evolution of the ALD growth surface.For example, [26] developed a 1D reactor model comprised of the species conservation equation to describe precursor concentration through the reactor vessel and a surface coverage equation, which relates adsorption to an experimentally-determined sticking coefficient based on time-dependent surface coverage and precursor concentration at the reactor outlet.The model was further used to study temperature effects on surface coverage [27] and the effects of secondary reactions on film thickness [28].A 1D plug-flow model is presented by [29], which also employs the species conservation equation to determine the spatially-and time-dependent pressure profile in the reactor starting with solid precursor evaporation, and uses a kinetic expression based on an estimated sticking coefficient.A simple plug flow model, based largely on empirical sticking and recombination probabilities, also was developed by [30] to study film conformality as a function of aspect ratio, but the details of the mathematical model were omitted.
A transient plug-flow model for a tubular ALD reactor was developed by [31], where the continuity equation is used to describe transport and the reactive sticking coefficient was based on earlier work of the author [32].Additionally, see [33] for a similar approach and a more general discussion of ALD reactor modeling.Furthermore, [34] describe a dynamic reactor model for alumina ALD processing that combines precursor flow modeling through the reactor and a ballistic flux model for precursor transport and reaction in substrate trenches.In that study, the reactive sticking coefficient was attributed to quantum-chemical simulations of the ALD surface reactions, but the connection between the two was not made clear.
Up to this point, the ALD transport/reaction models discussed have incorporated a reactive sticking coefficient or reaction probability, combining the quantitative effects of many reaction phenomena into a single parameter, which may be determined empirically or theoretically.Although the sticking coefficient is widely used in modeling work, the experimental and computational manner in which it is determined varies, such that absolute values are rarely reported in the literature [35].An alternative to the sticking coefficient approach is to couple transient plug-flow reactor dynamics with kinetic expressions containing rate constants derived from ab initio quantum-chemical calculations and the quantum-statistical theory of chemical reactions [36].Yet another approach is to use the absolute reaction rate theory and statistical thermodynamics to derive kinetic expressions without the use of adjustable parameters [11].Furthermore, [37,38] developed a rigorous 2D transport model, which is coupled to a heterogeneous surface reaction model based on estimated kinetic parameters from ex situ film thickness measurements.In this two-part paper series, the researchers examine film growth and thickness uniformity as a function of process operating parameters, reactor system design and gas flow distribution as a guide for future ALD process optimization.

The Reactor
The ALD reactor system considered in this work is based on a laboratory-scale research reactor currently under construction and to be used in evaluating ALD for a range of spacecraft-related thin-film applications [39].As illustrated in Figure 3, the reactor vessel consists of a stainless steel process tube surrounded by a bench-top tube furnace, containing the substrate(s) and a quartz crystal microbalance (QCM) for real-time mass deposition measurements.The reactor performance will be initially tested with the commonly used ALD precursors of trimethylaluminum (TMA) and water; both are contained in the liquid state in temperature-controlled laboratory bottles.Each of the precursors flows through a separate sequence of needle valves/orifices to control their flow rates and, then, through solenoid-activated control valves, CV1-5, to regulate the precursor dosages and the lengths of the purge periods.The details of the valve sequencing will be described in the Process Recipe section that follows.A mass flow controller (MFC) is used to regulate the purge gas (Argon) flow rate, and low-pressure manometers (denoted P1 and P2) respectively record gas pressure at the reactor outlet and the small gas chamber used to regulate the TMA dose.Residual gas analysis (RGA) is performed using a mass spectrometer; the primary vacuum pump is located downstream of the RGA, and a smaller pump is used to vent water vapor between water doses.Reactor dimensions and other reactor component specifications are listed in Table 1.
Table 1.Reactor dimensions and primary system component design parameters of the ALD reactor system shown in Figure 3. TMA, trimethylaluminum.

Parameter Description Value
A s substrate growth surface area The overall reaction involves trimethylaluminum Al(CH 3 ) 3 (precursor A) and water (precursor B) producing the amorphous Al 2 O 3 film and by-product methane gas; the CH 4 does not react with any of the other species in the deposition process and, so, can be considered inert.As described in the Introduction, each precursor is introduced sequentially, separated by a purge period to prevent gas-phase reactions: where the AB cycle repeats-potentially thousands of times-building the film one sub-monolayer at a time.After the initial nucleation transient following a change in the precursor system (e.g., when depositing a nanolaminate consisting of alternating thin-film materials), the deposition behavior during each AB cycle approaches a limit-cycle solution, the computation of which is the main focus of this simulation study.The exposure time periods (in seconds) for the AB exposures and purge periods are denoted as τ A , τ B , τ AP and τ BP .A nominal A-purge-B-purge process recipe will be described later in this manuscript.
Referring to Figure 3, control valves, CV2 and CV5, normally are open during all exposure and purge periods and, so, will not be discussed further.During the purge period prior to Exposure A, CV1 is closed to allow TMA to fill the ballast chamber; the TMA partial pressure in this chamber can potentially reach the vapor pressure of TMA in the supply bottle containing liquid-phase TMA.During Exposure A, CV1 is opened to allow TMA vapor to flow into the reactor, reducing the pressure in the ballast chamber.We note that a small flow of TMA will continue through the orifice/needle valve during Exposure A. At the end of Exposure A, CV1 is closed, and the pressure rebuilds in the TMA ballast chamber.
Regardless of the position of CV4, Ar purge gas flows continuously during all purge and exposure periods.During both purge periods and during Exposure A, CV3 is open between the water source and the water purge pump, as well as to CV4; however CV4 is closed in the direction of CV3, resulting in no water flow to the reactor.During the water dose, CV4 is switched to the all-open position, but CV3 is closed in the direction of the water purge pump, allowing the flow of Ar and water to the reactor.This configuration was designed to prevent condensation in the water delivery system and to improve the reproducibility of the water dose [40].

Precursor Characteristics
The ALD reactor system model development begins with the precursor thermophysical property and gas delivery system dynamics modeling.From the National Institute of Standards and Technology (NIST) WebBook [41], we find the Antoine's equation coefficients for TMA (between 337-400 K) and water (between 293-343 K) as: where the vapor pressure units are Torr and T is in K.The vapor pressures calculated by these relationships corresponding to the TMA and water sources are given in Table 2.As pointed out in [42], at lower temperatures and higher TMA partial pressures, the dimer of TMA (d-TMA) is favored over the monomer (m-TMA) in the gas phase.This is important in determining the true dose values for the precursor delivery system.The d-TMA dissociation extent as a function of temperature was studied in [43], where values for the degree of dissociation, α, were given as a function of temperature; the degree of dimer dissociation is defined as follows: where n o A is the total number of Al(CH 3 ) 3 molecules whether in monomer or dimer form.Using a least-squares procedure, we fit the following expression to the data contained in the cited source to find: where T again is in K.For a binary mixture of d-TMA and m-TMA, we can solve: for the mole fraction, y m , of the monomer with P o = 760 Torr.We note that the m-TMA mole fraction is related to the degree of dissociation, α [43], by y m = 2α/(1 + α) when considering a binary mixture.TMA bubbler/ballast chamber flow coefficient In a mixture containing additional species that do not participate in the dimerization reaction, defining n N as the total moles of species not participating in the d-TMA/m-TMA reaction, the gas-phase m-TMA and d-TMA mole fractions, y m and y d , can be written in terms of α, as follows: allows us to compute: Given the conditions inside the TMA source (T = 300 K, P = P vap T M A and φ = 0), the extremely small values of K d and α listed in Table 2 show that TMA leaves the bubbler essentially entirely as the dimer, d-TMA.This can also be observed in Figure 4, where the TMA vapor pressure and degree of dissociation are plotted.The objective of this model element is to predict the time-dependent TMA molar flow rate, ṁA in (t), as this precursor enters the reactor chamber.As seen in Figure 3, the precursor delivery system is designed to inject a TMA dose regulated by the size of a ballast chamber and the TMA pressure, P 2 , in this chamber prior to the opening of control valve, CV1.The dependence of gas molar flow rate through CV2 and the downstream orifice/needle valve would be a significant challenge to predict by a purely physically-based modeling approach.Therefore, under the condition that the gas flow is not in the choked-flow regime, we propose the model: where α source and α bc are the degrees of d-TMA dissociation in the TMA source and ballast chamber, respectively, computed using Equation ( 2) with φ = 0.The reactor and ballast chamber pressures are: where P rxr is the reactor pressure, approximately measured by manometer P1.The function, γ 1 , indicates the time-dependent position of CV1.We note that it is possible to measure the effective flow coefficients, C 1 and C 2 , experimentally using time-dependent measurements of P 2 .CV2 is always open under typical operating conditions, allowing for finite gas flow rate whether or not CV1 is open.We note that for this model, any back-diffusion of reactor Ar or other gas species is neglected.

Water Delivery System Model
Because of the potential for condensation of water in a ballast chamber for this precursor, the alternative design developed by [40] is used where water evaporating in the H 2 O source is vented using an auxiliary purge pump during the TMA dose and purge periods (a similar approach cannot be used for TMA, because of the expense of discarding unused precursor).This configuration gives rise to a relatively simple model for the water dose: The function, γ 3 , indicates the time-dependent position of CV3.Methods for measuring ṁB in are described in [40] and, so, can be used to identify the effective flow coefficient, C 3 .

Reactor Model
The definitions of the instantaneous consumption rates of the precursors, TMA, Γ A , and water, Γ B , and methane production, Γ C , per unit area of the growth surface, due to the deposition reactions, will be derived in Section 4.5 as Equations ( 21)- (23).Defining m rxr (t) as the total moles of gas-phase species in the ALD reactor, ṁin , the total reactor feed molar flow rate, and ṁout , the molar flow of residual gas out of the reactor, an overall reactor material balance gives: Vacuum pump capacity is rated by the volumetric pumping capacity per unit time, S t , given in Table 1.
Over the pump's operating range, this value is considered pressure-independent. Based on m rxr , the reactor instantaneous total pressure is computed using the ideal gas law Equation (4); so, the residual gas molar flow rate can be computed by: The total molar feed to the reactor is the sum of the precursor and inert flow rates, and so: where α rxr is the d-TMA dissociation degree under reactor pressure and composition conditions and is computed using Equation ( 2) and φ rxr .Under steady flow conditions with no precursor feed, such as at the end of a very long purge period, no surface reactions take place, and so, the base-line reactor pressure is defined using ṁout = ṁin = ṁI in and: in which for a nominal argon flow, ṁI in , and assuming that the residual gas has cooled to ambient temperature, T amb , by the time it reaches the pump, gives P rxr,base , listed in Table 2. Total reactor molar capacity (m rxr ), residence time (τ res ) and other parameters evaluated at the nominal conditions also are given in Table 2.
Defining y A as the gas-phase mole fraction of m-TMA plus d-TMA and y B , y C and y I as, respectively, the gas-phase water, methane and argon mole fractions, the reactor dynamic material balances can be written as follows: subject to initial conditions, y A (0), y B (0), y C (0), y I (0), m rxr (0), the time-varying state of the growth surface and the total molar inlet flow ṁin , given by Equation (6).

Surface Reaction Mechanisms and Kinetics
During Exposure A, the gas-phase m-TMA molecules (A) adsorb onto surface hydroxyl groups (X) or surface oxygen bridges (O), the latter corresponding to alumina film oxygen atoms located on the growth surface.In the first case, a Lewis acid/base reaction results in a chemisorbed surface adduct species (AX) comprised of a TMA and hydroxyl group, a mechanism examined by the DFT studies of [44].The second case results in a dissociative adsorption reaction, where TMA adsorbs onto an oxygen bridge and breaks an Al-C bond, leaving three Me groups on the growth surface [45] (we will not consider this reaction in this study).Subsequent to either of these adsorption steps, H migration and reaction with surface Me groups releases CH 4 , which desorbs immediately [10], effectively resulting in an irreversible reaction and depopulating the surface of hydroxyl groups.The reactions taking place during the water exposure follow structurally similar reaction pathways and are described in more detail in [11].Further details also are given in the excellent review by [13] of atomic-scale ALD reaction mechanisms and the quantum-chemical computations used to uncover the reaction mechanisms.

The Surface State
To characterize the growth surface, we denote [X], [O], [S], [AX], [BS] and [M e] as the hydroxyl, oxygen, aluminum, TMA-OH adduct, water-Al adduct and methyl group surface concentrations (species/nm 2 ), respectively.Despite the amorphous nature of the alumina film, we represent the instantaneous state of the surface in the manner shown in Figure 5, where the X,O,S checkerboard pattern corresponds to a grid of 0.295 nm × 0.295 nm squares and the Me group radius of 0.2 nm [11].How the maximum surface density values, [ X], [ Ô], [ Ŝ], [ M e], were computed also is described in the cited reference.Using these definitions, we can now define the surface coverage fractions.
Figure 5.A snapshot of a 35 nm 2 portion of the ALD growth surface corresponding to θ O = 0.6 and θ M e = 0.5.

Surface Reaction Equilibria and Reaction Rates
The Exposure A reactions between TMA and a surface hydroxyl group follow a sequence of adsorption, adduct formation, transition state and irreversible final reaction steps; the reaction mechanism postulated by [44] can be written as reaction R0:

X+A
adsorbed adduct AX transition state AX ‡ products where 0 ≤ 0 and ‡ 0 ≥ 0 are the adsorption and single irreversible reaction activation energies, respectively, and K 0 and K ‡ 0 are equilibrium coefficients.In this reaction sequence, the rate-determining step is the final irreversible reaction that adds Al(CH 3 ) 2 to the growth surface and produces one methane molecule that desorbs to the reactor gas phase.The two equilibrium reactions define the surface concentrations of the adduct, AX, and the transition state, AX ‡ , between the adduct and the final reaction products.
A rate expression for v 0 was developed in [11]; for this study, we consider the alternative reaction mechanism of [10], written as the following reaction sequences: R0 and R1 begin with the same TMA adsorption process: with A as the gas-phase TMA species and X the surface sites to which the precursor can be adsorbed to form surface adduct species, AX: assuming ideal gas behavior for A. 1 ≤ 0 is the potential energy change associated with TMA adsorption and adduct formation.The reaction energies for each reaction and the sources of those values are given in Table 3. Z AX , Z A and Z X are the adduct, gas-phase TMA and surface OH group partition functions, respectively.Solving for [AX] and considering all thermodynamic quantities to be on a per-molecule basis, we recover the Langmuir isotherm: Moving on to the adduct/transition-state equilibrium, AX + X AXX ‡ of R1, we observe that the primary difference between R0 and R1 is that the latter forms a transition state by reacting with a neighboring surface H rather than the hydrogen associated with the original OH adsorption site.The equilibrium relationship is written as: with ‡ 1 = ∆E P ≥ 0 the activation energy or change in potential energy associated with going from the adduct species to the transition state.Z AXX ‡ is the partition function for the transition state, and the [ X] is included as a means of approximating an infinitely slow surface diffusing X [46].We note that AXX ‡ has one fewer vibrational modes relative to the maximum number possible for this species, because of its role in the rate-limiting reaction [47], and this is reflected in the definition of Z AXX ‡ .Combining the adsorption and reaction steps, the final R1 reaction sequence rate expression is found to be: R4 activation energy 0.91 [44] Each of the two surface Me groups left after R1 can undergo subsequent reactions with surface OH groups by the mechanism depicted as R2.One can immediately write: giving the reaction rate: We note that surface oxygen sites, O, are produced by every X consumed in R1 and R2.The partition functions for the equilibrium coefficients are as follows: where s = A, X, AX, AXX ‡ , Me or MeX ‡ and Z vib s , Z rot s and Z trans s are the vibrational, rotational and translational components of each partition function, respectively.These partition function components are described in detail in [11].Z A has units of m −3 , and all other partition functions are dimensionless.

Water Reactions
Following [44], we write the water-exposure reactions as: Because reactions R3 and R4 essentially follow the same sequence of the equilibrium-limited adsorption and irreversible reaction steps of R0 and R1, we can immediately write: with i = 3, 4 and partition functions that are likewise similar to those of the TMA reactions.The difference between R3 and R4 is related to the surface OH concentration: R4 is more likely to take place with increasing X concentration relative to R3.We note, however, that the Al of reactions R3 and R4 has a coordination number of three, while R1 produces Al with a coordination number of four; as discussed in [10], the latter is energetically more favorable relative to the three-coordinated Al.Reconciling these surface reaction details will be the subject of a follow-up study.

Surface State Dynamics
Surface Me are produced by R1 and are consumed by R2, R3 and R4.Likewise, surface X are converted to O by R1 and R2, while X is created by R3 and R4.The essential self-limiting behavior of ALD processes results from surface Me saturation ([M e] → [ M e]) shutting down reaction R1.Combining these effects, we can write the surface species balances for Me and X as: subject to initial conditions, θ M e (t = 0) and θ X (t = 0), at the start of each exposure and purge period.

Gas-Phase Species Flux at the Growth Surface
We see that the surface phase rate equations are coupled to the gas phase through the precursor partial pressure, P A and P B , in Equations ( 15) and (18), respectively (Note that at this time, it is unclear whether d-TMA or only m-TMA can participate in the adsorption reactions, and so, the total pressure, P A , is used).To compute the rate of gas-phase depletion of the precursors, due to the surface reactions and the rate of production of the methane by-product, we need the rates of consumption (positive quantities) of TMA, water and methane, denoted as Γ A , Γ B and Γ C , respectively: molecules/(nm 2 s) ( 21)

Limit-Cycle Computations
At this point, we write the complete model as: and collect the differential equation right-hand sides, model variables and process recipe parameters in the following table: 3), ( 7)-( 10), ( 19), ( 20 where y I = 1 − y A − y B − y C and the length of the full process cycle is: All of our simulations are implemented in the Python programming language, making extensive use of the PyLab (www.scipy.org/PyLab)and Numdifftools (pypi.python.org/pypi/Numdifftools)modules.Therefore, any computationally specific discussions that follow will be in the context of a Python implementation.

Time Discretization
Our solution strategy for the dynamic ALD process is to only consider the limit-cycle solutions that describe steady (but periodic) operation of the reactor system.Nucleation and other events leading to transients spanning multiple exposure cycles are not considered in this work.Computation of limit-cycle solutions naturally lend themselves to a two-step procedure, where the modeling Equation ( 24) is first discretized in time over each exposure and purge period using global collocation over each of the four periods (τ A , τ AP , τ B , τ BP ), enforcing continuity between each interval, effectively resulting in periodic boundary conditions between the end of the BP purge and the start of the next A dose.The resulting nonlinear equations, then, are solved using the Newton-Raphson method.
To implement the collocation procedure, we must first define the format of the Python array used to represent the time-discretized vector of state variables, ξ ∈ R n , as Ξ ∈ R mn .For reasons advantageous to computing the Jacobian array elements in the Newton procedure, we define the Python list, Ξ, as a list of the process variables, where Ξ(i, j) is state j at point i in time; defined in this manner, Ξ(i) is a list representing a snapshot of the states at a specific time.Given this format for the discretized states, we can write the discretized form of Equation ( 24) as: where the Ξ list is flattened to the shape appropriate for matrix multiplication using the Python ravel method, and Â is defined below.With A n×n corresponding to the standard Lobatto first-order differentiation array (computed using either finite differences or using polynomial collocation techniques), the discretization array suitable for vectors of discretized states, Â, is created from diagonal m × m arrays from elements of A: Note that the identity array, I 1,1 , has dimensions n × n and is used to satisfy the initial conditions.In this study, each of our discretized intervals uses n collocation points (including both endpoints).Writing the vector of discrete points in time over each exposure and purge period: where the off-diagonal blocks, C, are used for continuity across the spline point, and the off-diagonal block P is used to enforce periodicity.

Newton-Raphson Procedure
With the discretization complete, we write the Newton-Raphson procedure in terms of the residual, r, update, u, and refined solution estimate, Ξ, at iteration ν: While perfectly standard, we present the Newton-Raphson procedure to point out the structure of the Jacobian array.Numerical approximation of the full Jacobian array does not take advantage of its structure: This limits the extent to which a finite-difference procedure must be applied to compute the block-diagonal Jacobian elements corresponding to the relatively complicated, nonlinear terms in the rate expressions, precursor thermodynamics descriptions and reactor material balances.The Jacobian array, J, then, is constructed in a block-diagonal manner, calling the Python function, numdifftools.Jacobian, for each (collocation) point in time to define the block-diagonal elements of Equation ( 27).

Results
A limit-cycle solution is presented in Figure 6 corresponding to a base-case design, both in terms of the reactor component specifications and the process recipe.The nominal design parameter values have been listed in Tables 1-3.Three sets of plots are presented in Figure 6 in which subplots illustrating the TMA ballast chamber, reactor gas phase and growth surface composition dynamics are shown.One can observe that all states conform to periodic boundary conditions over the processing cycle; the markers indicate the locations of the temporal collocation points and the shaded rectangles, the collocation interval endpoints.

TMA Ballast Chamber Dynamics
The full process limit-cycle is shown to begin with Exposure A, where the valve between the ballast chamber and reactor is opened; the valve between the TMA source and ballast chamber (CV2) always is open and has a fixed flow coefficient through the entire process cycle.During dose A, the TMA flows from the ballast chamber to the reactor chamber, reducing the pressure of the former.At the end of Exposure A, the ballast chamber/reactor valve, CV1, is closed, allowing the TMA pressure to rebuild during the subsequent purge and Exposure B periods.
As seen in Figure 6, top, the ballast chamber pressure behaves exactly as expected, with a rapid initial drop in pressure, due to the small volume, V bc .However, what is interesting to observe is the degree of d-TMA dissociation α bc in the bottom plot of Figure 6: except for a slight upward deviation during the ballast chamber depressurization, α bc ≈ 0, indicating that virtually all of the TMA is in dimer form while in the TMA ballast chamber.

Reactor-Scale Dynamics
As seen in the center plot of Figure 6, during Exposure A, the total reactor pressure, TMA partial pressure and methane partial pressure all increase as expected, while the Ar carrier gas partial pressure remains constant.During the subsequent purge period, the total pressure relaxes to the base-line pressure.It is now interesting to observe that the TMA monomer fraction in the reactor is essentially unity (α rxr ≈ 1), indicating that the d-TMA dissociates as it enters the lower-pressure, higher-temperature reactor chamber.We note that the energy required to heat the incoming reactant and inert gases is negligibly small compared to the rates of radiative heat transfer in the reactor, and so, we do not explicitly model the thermal dynamics of the gases as they are heated to T rxr from T amb or T bc .
During Exposure B, we observe a much more linear increase in total pressure, because water is fed to the reactor from a water vapor source held at constant pressure with P vap water P rxr,base .Again, the system relaxes to the base-line pressure, P rxr,base , after valves CV3 and CV4 are switched to their purge positions.

Growth Surface Dynamics
At the start of dose A in the limit-cycle solution, θ M e ≈ 0.08 and θ X ≈ 0.25.As TMA enters the reactor, a small fraction rapidly reacts with the surface OH, leading to a reduction in θ X and the increase in θ M e shown in Figure 6.With sufficient TMA dosage levels, θ M e → 1 very rapidly, indicating the aggressiveness of reaction R1.As the growth surface saturates with Me groups, the rate of R1 slows, and R2 becomes the rate-controlling step.We observe that R2 continues throughout the purge period, reducing both the surface Me and OH ligand density.
When H 2 O is introduced during Exposure B, θ M e rapidly drops as the surface Me groups are replaced with OH in reactions R3 and R4.As the water partial pressure rapidly drops during the post-B purge, all reactions, except R2, come to a stop, leaving a nearly unchanging growth surface for much of the second purge period.The full length of the purge period is required, however, to prevent remaining gas-phase H 2 O from reacting with surface Me once the dose-A period resumes.

Mapping the gpc Behavior
As described earlier, a distinguishing feature of ALD processes is the self-limiting nature of the surface reactions, leading to a fixed rate of growth-per-cycle (gpc) during steady, but cyclic, reactor operation.The total number of Al and O atoms deposited per unit area over one deposition cycle are denoted as N Al and N O atoms/nm 2 , respectively; these values are computed by integrating the consumption rates of both precursors, Equations ( 21) and (22), over the limit-cycle: using the quadrature weights of the collocation procedure.The relative ratio of Al to O, thus, can be determined; for a ratio of Al/O = 2/3: where ∆ z is determined from the density of amorphous Al 2 O 3 [11].Alternatively, if η 1 is the extent of reaction R1, the only reaction involving gas-phase TMA, we can write: where the subscripts, o and f , denote surface concentrations at the start of dose A and the end of post-A purge, respectively.We note this relationship holds only when no reactions take place under CVD conditions (where both gas-phase precursors are found in the reactor, resulting in the possibility of reactions R1-R4 all taking place simultaneously).Under these idealized ALD reactor operating conditions and under fully saturating conditions, resulting in Equation ( 30), reducing to the maximum gpc possible for an idealized ALD process, a value denoted as GP C: We note the value of GP C correlated well with previous theoretical and experimental studies of the alumina ALD process [48].
The gpc = 0.806 Å/cycle of our base-case design corresponding to the limit-cycle solution of Figure 6 is considerably less than the maximum indicated by Equation (31) for idealized, saturating ALD conditions.We now examine two modes of increasing the dose of each precursor.During Exposure B, the water dose is regulated by the timing of the gas delivery system valves, CV3 and CV4.With a base-case design of τ B = 0.1, we observe in Figure 7, left, that gpc → 0 as τ B → 0, while keeping V bc fixed, exactly as expected.We note that τ B = 0 actually corresponds to a bifurcation point, where the branch containing the limit-cycle solution shown in Figure 6 meets a trivial solution characterized by θ X = θ M e = 0 for all points in time.The physical meaning of the multiple solutions will be explored in follow-up work.As τ B is increased from the nominal operating conditions, the rate of gpc increase lessens; the CVD conditions and the slower surface reactions contribute to the gradual increase in gpc with no self-limiting regime to be found under the selected set of operating conditions of the plot.The base-case ratio of TMA ballast chamber/reactor chamber volume is 0.02%, and one expects that increasing this ratio will result in an increased TMA dose in the reactor vessel.Keeping the H 2 O dose fixed, the overall results are seen in Figure 7, right.As V bc → 0, gpc → 0.78 Å/cycle-not zero-because τ A is nonzero, and the TMA bleeding through valve CV2 always results in a nonzero TMA dose.We observe that gpc grows with V bc , but again, while the rate of gpc increase declines with ballast chamber volume, a plateau indicating saturating growth is not observed under the operating conditions of the plot.

The V bc -τ B Plane
Because the gpc values of Figure 7 neither reach a limiting value nor the predicted theoretical maximum, gpc → GP C = 1.231Å/cycle, we present the gpc as the contour plot in the V bc -τ B plane in Figure 8.In this figure, the base-case design corresponds to the lower-left corner of the plot.In Figure 8, we mark gpc = GP C = 1.231Å/cycle by the white contour line found over this range of τ B and V bc values.The large region above and to the right of the curve corresponds to large doses of both TMA and water, resulting in growth exceeding what would be possible for a pure, surface reaction-limited ALD process.Examining the limit-cycle solution in this region reveals that because the amount of TMA and water is so large, a significant amount remains after each purge period, resulting in each of the precursors being found in the gas phase at the start of the other dose period.During Exposure A, excess water in the gas phase reacts with surface Me deposited by TMA, increasing the ability of TMA to adsorb and react during dose A. Likewise, excess TMA present at the start of dose B generates surface Me, which subsequently react with gas-phase H 2 O, adding to the overall consumption of that reactant.
These unintentional reactions are characterized as being under CVD conditions instead of true ALD reactions.Because CVD conditions are generally undesirable in ALD processes, due to the poor film quality and reduced produced by the resulting reactions, we mark the curve corresponding to gpc = GP C to indicate an approximate lower limit, where reactions under CVD conditions are significant.Thus, ALD reactor operation should be limited to the region below and to the left of this curve.The practical upshot of this computation is immediate: we see there is little incentive for process designs where V bc /V rxr > 0.25% and that, generally, τ B > 0.2 s (given that all other parameters are fixed, of course).The rationale for these limits is further clarified by considering the economics of this ALD process: given the relative value of TMA to water, plus the costs of disposing unused TMA downstream of the reactor, a simple yet reasonable optimization objective would be to minimize TMA use alone.To quantify TMA use, we integrate the TMA flow rate through CV1 using the ṁA in term of Equation ( 6 These curves are shown in black in Figure 8, where the values indicated correspond to moles of TMA/cycle fed to the reactor.As expected, the contours show a reduction of TMA use as V bc is decreased, but we note that for large dose volumes, both τ B and V bc affect m A cycle , due to the increased time for regenerating the TMA pressure in the ballast chamber during dose B.

Conclusions
In this paper, a laboratory-scale atomic layer deposition (ALD) reactor system model was developed by integrating components describing the precursor thermophysical properties, reactor-scale gas-phase dynamics and surface reaction kinetics derived from absolute reaction rate theory.ALD reactor operation was limited to steady cyclic operation with limit-cycle solutions computed using a collocation discretization scheme in time.We demonstrated that a key advantage to the fixed-point approach was the resulting unambiguous definition of growth-per-cycle (gpc), making possible parametric studies of film growth rates to that expected for ideal ALD.
The utility of the resulting ALD system simulator was demonstrated by the strong interactions found between different reactor, reaction and process recipe elements, interactions that would be otherwise difficult to predict without such simulators.In particular, we demonstrated that surface reactions normally associated with one specific precursor exposure can take place during the purge or even during exposures to the other of the two precursors.The ability to model the interaction between dose and purge periods was critical to uncovering the surface reactions occurring under CVD conditions and identifying processing regimes where these reactions are most likely to take place.
The ability to predict both gpc and cyclic precursor feed rates for a real reactor and gas delivery system made it possible to generate simple, but physically meaningful, design rules for adjusting the precursor doses to minimize TMA consumption and undesirable CVD conditions, while maintaining a high value of gpc necessary for acceptable reactor throughput.One of the most important contributions of reactor system-level models for thin-film processes is the ability to use dynamic models to accurately characterize the time-dependent composition of the reactant gases to which the growth surface is exposed.Because direct experimental measurement of gas-phase characteristics local to the growth surface are challenging at best, physically-based reactor models are needed to interpret measured gpc levels.The utility of models of this form extend to other ALD process chemistries, other gas delivery systems and more complex (e.g., multi-wafer) reactor systems.The extension of our modeling work to the tubular geometry of our ALD system is underway.Likewise, the use of our reactor models in more sophisticated dynamic optimization studies, as well as controller development to fully decouple binary precursor doses are planned.

Figure 1 .
Figure 1.An idealized view of the atomic layer deposition (ALD) process cycle.

Figure 2 .
Figure 2. Elements of a complete ALD reactor system mathematical model.

Figure 3 .
Figure 3.A schematic diagram of the laboratory-scale ALD reactor system to be modeled.

Figure 4 .
Figure 4. TMA vapor pressure and degree of dissociation, α, as a function of temperature; data presented in TableIIof[43] are shown as filled red circles.

Figure 7 .
Figure 7. Growth per cycle (gpc) as a function of τ B (left) and V bc /V rxr (right).

Figure 8 .
Figure 8. Growth per cycle (gpc) as a function of τ B and V bc /V rxr .The black curves correspond to moles of TMA/cycle fed to the reactor.

Table 3 .
Reaction rate energetics information.