A model-based methodology for the analysis and design of atomic layer deposition processes—Part II:

Chemical Engineering Science 94 (2013) 316–329

Contents lists available at SciVerse ScienceDirect

Chemical Engineering Science journal homepage: www.elsevier.com/locate/ces

A model-based methodology for the analysis and design of atomic layer deposition processes—Part II: Experimental validation and mechanistic analysis b ¨ ¨ a A. Holmqvist a,n, T. Torndahl , S. Stenstrom a b

Department of Chemical Engineering, Lund University, P.O. Box 124, SE-221 00 Lund, Sweden ˚ Angtr¨ om Solar Center, Solid State Electronics, Uppsala University, P.O. Box 534, SE-751 21 Uppsala, Sweden

H I G H L I G H T S c c c c c

Methodology for estimating kinetic parameters involved in the ALD gas–surface reactions. Assessment of the statistical reliability of the parameter estimates and the model response. Continuous flow ALD reactor model was validated experimentally and showed high predictability. Mechanistic dependence of process operating parameters on ALD film growth profiles were assessed. Application of a heuristic optimization algorithm implemented in a cluster environment.

a r t i c l e i n f o

a b s t r a c t

Article history: Received 20 February 2012 Accepted 28 June 2012 Available online 11 July 2012

This paper demonstrates the experimental validation and mechanistic analysis of the continuous crossflow atomic layer deposition (ALD) reactor model developed in the first article of this series (Holmqvist et al., in press). A general nonlinear parameter estimation problem was formulated to identify the kinetic parameters involved in the developed ALD gas–surface reaction mechanism, governing ZnO film growth, from ex situ film thickness measurements. The presented methodology for comprehensive model assessment considers the statistical uncertainty of least-squares estimates and its ultimate impact on the model predicted response. Joint inference regions were determined to assess the significance of parameter estimates and results indicate that all estimates involved in the precursor half-reactions were adequately determined. The reparameterization of the Arrhenius equation effectively decreased the characteristically high correlations between Arrhenius parameters, leading to improvement in precision of individual parameter estimates. Model predictions of the spatially dependent film thickness profile with narrow confidence band were in good agreement with both calibration and validation experimental data, respectively, under a wide range of operating conditions. The subsequent extensive theoretical analysis exhibits that the experimentally validated model successfully reproduces the detailed process dynamics revealed by in situ quartz crystal microbalance and quadrupole mass spectroscopy diagnostics, and thereby provides a supplementary analysis tool. Finally, the univariate sensitivity analysis revealed the mechanistic dependence of all the measured process operating parameters on the spatially dependent film thickness profile, resolved at the level of a single pulse sequence. Hence, the presented model-based framework serves as a means to guide future research efforts in the field of ALD process optimization. & 2012 Elsevier Ltd. All rights reserved.

Keywords: Atomic layer deposition Experimental model validation Parameter identification Optimisation Dynamic simulation Kinetics

1. Introduction Knowledge of the microscopic physicochemical growth process and the underlying macroscopic mass transport is a fundamental

n

Corresponding author. Tel.: þ46 46 222 8301; fax: þ46 46 222 4526. E-mail address: [email protected] (A. Holmqvist).

0009-2509/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ces.2012.06.063

requirement for the understanding and modeling of atomic layer deposition (ALD), see part I of this article series (Holmqvist et al., in press or Elliott, 2007). According to the International Technology Roadmap for Semiconductors (ITRS) (Semiconductor Industry Association, 2011), process models have the potential to yield insight that will enable improved efficiency and throughput in manufacturing and in the design and evaluation of reactors. Due to the lack of reliable fundamental physical and chemical data for new

A. Holmqvist et al. / Chemical Engineering Science 94 (2013) 316–329

precursors, and the introduction of new process configurations and regimes, calibration and experimental validation of process models have been identified as an important step in the field. There are two main ways of providing additional physical insight into the nature of the reaction mechanisms concerned and determining reliable estimates of the kinetic parameters. (i) Ab initio quantum-chemical and quantum-kinetic calculations (Deminsky et al., 2004; Xu and Ye, 2010), and (ii) model calibration and experimental validation (Lim et al., 2000; Park et al., 2000). However, the main objective of the published work was to determine the principal features of the kinetic mechanism, and detailed fluid dynamics was neglected at that stage. It was emphasized by Puurunen and Vandervorst (2004) that model calibration utilizing experimental saturation curves determined from flow-type reactors, which are dominated by mass transport (Aarik and Siimon, 1994; Siimon and Aarik, 1995; Ylilammi, 1995), are generally not suitable for extracting kinetic parameters governing ALD growth without a reactor scale model. For this reason, the novel mathematical model of the continuous cross-flow reactor system F-120 by ASM Microchemistry Ltd. (Suntola, 1992) presented in the first article of this series (Holmqvist et al., in press) was derived to mechanistically describe the complex interactions between the gas-phase fluid dynamics, the mass transport of individual species, and the ALD heterogeneous ¨ 2002). Thereby, surface reaction mechanism (Ritala and Leskela, the model incorporates all relevant physicochemical phenomena governing ALD growth required to assist in parameter estimation. Furthermore, in contrast to the mathematical models proposed in the aforementioned cited studies this model, comprises a system of nonlinear partial differential equations (PDEs), provides spatial resolution of the substrate film thickness profile (Siimon and Aarik, 1995, 1997). Given the novel application, the detailed experimental investigation on growth of ZnO films from Zn(C2H5)2 and H2O reported in Holmqvist et al. (in press) was performed in order to assess the validity of the derived mathematical model on a discrete set of sampling positions on the substrate, see e.g. Cleveland et al. (2012) and Henn-Lecordier et al. (2011). The overall objective of this study was to experimentally validate and demonstrate the generally applicable mechanistic model for the analysis and design of ALD processes developed in the first article of this series (Holmqvist et al., in press). This paper focuses around three main points. (i) To formulate and solve a nonlinear parameter estimation problem using ex situ experimental data and the modeling environment presented in Holmqvist et al. (in press), and (ii) to assess the statistical reliability of the parameter estimates and the predicted model response, and (iii) to perform an extensive analysis of the experimentally validated ALD process model dynamics, and to mechanistically determine the sensitivity of the equipment-scale process operating parameters on film growth. Although assessed in several experimental studies, see e.g. Aarik et al. (2006), Jur and Parsons (2011), and Lei et al. (2006), the mechanistic details of the process operating parameters have not yet been thoroughly studied theoretically. The reminder of this paper is organized as follows. The nonlinear parameter estimation problem is formulated in Section 2, while the subsequent section describes the statistical model validation methodology. Section 4 outlines the optimization environment. Section 5 presents the results from the validation methodology followed by an extensive mechanistic analysis of the ALD process model. Finally, concluding remarks are drawn in Section 6.

2. Nonlinear parameter estimation The inverse method used in this study utilizes the fact that the height of the film contains information on the elementary surface

317

reactions and can thus be used to estimate the kinetic parameters (Holmqvist et al., in press):

b ¼ ½L,kT ref,i ,i ,Ei y

8i A f1; 2, . . . ,4g

ð1Þ

Nb

where b A R is the Nb -dimensional parameter space. However, identifying the restricted set of kinetic parameters, b, in an ALD process model represents a non-trivial task. The two main reasons are: (i) There is a complex interdependence between the sequential and parallel elementary surface reactions (Pinto et al., 2011) involved in ALD, such that the reactivity in one half-cycle is influenced by that in the half-cycle preceding it (Kuse et al., ¨ 2002). Consequently, the kinetics of 2003; Ritala and Leskela, each reaction step must be balanced with respect to the other. (ii) The intrinsic mathematical structure of the nonlinear Arrhenius equation introduces high correlation between the frequency factor, Ai, and the activation energy, Ei (Buzzi-Ferraris and Manenti, 2009; He´berger et al., 1987; Nagy and Tura´nyi, 2011; Schwaab and Pinto, 2007, 2008; Schwaab et al., 2008b; Varga et al., 2011). For this purpose, a carefully chosen experimental design have been recognized as an important element in comprehending the coupled nature of the ALD growth kinetics and to obtain adequately determined parameters (Park et al., 2000). The aim of the experimental design is to yield the most possible information, in a statistical sense, for use in parameter estimation and model validation (Franceschini and Macchietto, 2008). Further, given a set of discriminating experiments (Alberton et al., 2010, 2011), it is important to utilize the experimental information in an optimal preserving manner. As emphasized by Schwaab et al. (2008a), the error structure of the kinetic constants in nonlinear models may be non-convex, open and constitute of discontinuous regions. For this reason, it was recommended that the entire calibration parameter vector, b, should be calibrated simultaneously using ^ presented in the entire experimental calibration data set, y, Holmqvist et al. (in press). 2.1. Formulation of the nonlinear parameter estimation problem The nonlinear parameter estimation problem is formulated as a general nonlinear optimization problem: min bAR

Nb

FðbÞ

  dx ,x,u, b,w,t , 0¼F dt y ¼ gðx,u, b,wÞ bmin r b r bmax

s:t:

with

b ¼ ½L,kT ref,i ,i ,Ei y y

xE ¼ ½p,v, oa , yk ,ms  y

yE ¼ ½hs ðzÞ9f^ ,UF y 

8i A f1; 2, . . . ,4g

8a, k   ^f ¼ 1 , 3 , 5  9G 9 sub 6 6 6

uE ¼ ½Dt ZnðC2 H5 Þ2 , Dt H2 O ,Ty wE ¼ ½Dt N , Q_ , Q_ , Q_ 2

N2

xðt 0 Þ ¼ x0

ZnðC2 H5 Þ2

y H2 O ,/pS9Gout 

ð2Þ

where b is subject to lower and upper bounds acting as inequality constraints and estimated by minimizing an objective func^ and tion FðbÞ, based on the deviations between the observed, y, predicted system response, y. By using the method of lines (Davis, 1984; Schiesser, 1991) and finite element method (Zienkiewicz s and Taylor, 2000a,b) framework in COMSOL Multiphysics (COMSOL AB, 2008), the associated system of nonlinear differential algebraic equations (DAEs), in which x, u and w represent dependent states, free design variables, and algebraic variables, is

318

A. Holmqvist et al. / Chemical Engineering Science 94 (2013) 316–329

denoted by F in Eq. (2). It is noteworthy that response function g, which transforms and selects those state variables that are experimentally measured, is a nonlinear function of b for all states x. The weighted sum of squared residuals is used to quantify the estimation, and is defined as X ^ FðbÞ ¼ ½y^ E yE ðxE ,uE , b,wE Þy W1 ð3Þ E ½y E yE ðxE ,uE , b,wE Þ E¼1

where W is the square weight matrix corresponding to the different weights used to normalize the contribution of each particular observation E within the N E -dimensional sample space. The details underlying the statistical assumptions and their implications, which lead to the use of the least-squares estimates encompass several different aspects of the regression model, are presented elsewhere (Bates and Watts, 1988).

3.2. Assessment of the accuracy and reliability of the predicted model response The net effect of the statistical assumptions that are inherently made when applying the least-squares estimates and the propagation of parametric uncertainty is to produce model predictions with an uncertainty bound (Draper, 1995). The approximate 100ð1aÞ% confidence bands for the response function (Bates and Watts, 1988) are given by the collection of all curves yðx,u, b,wÞ belonging to the set: n yðx,u, b,wÞ : 9fyðx,u, b,wÞgfyðx,u, b^ ,wÞg9 r ½Jðb^ Þy Rðb^ ÞJðb^ Þ1=2 ðNb F ðNb , j;aÞ Þ1=2

o

ð7Þ

where Jðb^ Þ ¼ @yðx,u, b,wÞð@bÞ1 9b ¼ b^ denotes the parameter Jacobian matrix (Chappell and Gunn, 1998). 3. Statistical model validation via nonlinear parameter estimation In accordance with the methodology employed by Pinto et al. (2011), the statistical model validation methodology developed for this study is subdivided into two main parts: (i) parameter identifiability and the uncertainty of the estimates and (ii) goodness-of-fit of the model response and the underlying assumptions regarding the disturbances. 3.1. Assessment of the accuracy and reliability of the parameter estimates Parameter precision is equivalent in mathematical terms to the size of the inference regions of the parameter estimates (Hamilton et al., 1982). Based on the assumptions underlying least squares estimation, see e.g. Bates and Watts (1988), secondorder statistics (Witkowski and Allen, 1993) were explored for the construction of the confidence regions and assessing the uncertainty of b^ in the parameter space b A RNb . The approximate 100ð1aÞ% joint confidence region, obtained from a second-order Taylor expansion of the objective function, Fðb^ Þ, around the minimum point estimate (Schwaab et al., 2008a), is defined by fb : ½bb^ y Rðb^ Þ1 ½bb^  r Nb F ðNb , j;aÞ g

ð4Þ

where F ðNb , j;aÞ is the upper a quantile for Fisher’s F-distribution with j ¼ NE Nb degrees of freedom and Rðb^ Þ is the ðNb  N b Þ parameter variance–covariance matrix (Asprey and Macchietto, 2000). This approach also suggests the 100ð1aÞ% marginal confidence interval for bı (Bard, 1974; Draper and Guttman, 1995), ıA f1; 2, . . . ,Nb g: fb : 9bı b^ ı 9r Sıı ðb^ Þ1=2 t ðj;1a=2Þ g

ð5Þ

where t ðj;1a=2Þ is the upper a=2 quantile for Student’s t-distribution and Sıı ðb^ Þ is the ðı,ıÞ th element of Rðb^ Þ. Nevertheless, as emphasized by Bates and Watts (1988), Watts (1994) and Witkowski and Allen (1993), this linearization approach is only approximately applicable to nonlinear models, and in the case of high nonlinearity it can be quite inaccurate. For this purpose, the likelihood 100ð1aÞ% joint confidence region (Beale, 1960) is defined for all values of b such that fb : FðbÞFðb^ Þ rs2 Nb F ðNb , j;aÞ g

ð6Þ

where s2 ¼ Fðb^ Þj1 is the variance estimate of s2 . In contrast to Eq. (4), the confidence regions obtained with this method can be disjointed and unbounded, depending on the complexity of the parameter nonlinearity.

4. Optimization environment The minimization of the objective function, FðbÞ, in nonlinear parameter estimation problems is associated with complex numerical problems (Buzzi-Ferraris and Manenti, 2009; Moles et al., 2003). Therefore, the following measures must be taken into consideration in the formulation of the method: (i) the size of the parameter space, (ii) the multimodal nature of the objective function, (iii) the continuity of the objective function, and (iv) the sensitivity of the objective function to each of the model parameters (Hibbert, 1993). One way to overcome the non-convexity of the optimization problem is to consider a method of heuristic optimization, such as evolution strategies and genetic algorithms (Price et al., 2005). These algorithms overcome the starting-point problem resulting from the multi-modal nature of the objective function, by sampling the objective function at multiple, randomly chosen initial points where the present parameter bounds define the domain, and rely on a large number of objective function evaluations and a random search character to assure a high probability of finding the global optimum. These characteristics lead to prohibitively long computation times for the solution of this kind of problem, making nonsequential solution approaches necessary. According to the ITRS initiative, time-critical simulation steps or algorithms should support parallelization on high-performance computation cluss ters. Hence, the forward model simulator, COMSOL Multiphysics (COMSOL AB, 2008) integrated with MATLABs (The MathWorks Inc, 2010a) was extended for implementation in a cluster environment, following the master-slave paradigm (Coello et al., 2007). The model-based methodology supports the working procedure as follows: First, the parameter space, RNb , is sampled with Latin hypercube sampling (LHS) (McKay et al., 1979), using lhsdesign in TM the Statistics Toolbox for MATLABs (The MathWorks Inc, 2010c) to determine a feasible input for the heuristic optimization method. Subsequently, the evolution strategy differential evolution (DE) DE/ rand-to-best/2/bin (Storn and Price, 1997; Price, 1999; Price et al., 2005) was used to find the optimal least-squares estimates. Applying the concept of distributed computation resulted in a scalable and affordable cluster solution for the use of both the LHS and DE algorithms. Nonetheless, gradient methods outperform direct search methods both in reliability and speed of convergence (Bard, 1974). Additionally, DE possesses no strict convergence criteria. Therefore, in order to ensure reliability and to further enhance convergence of the output optimal calibration parameter vector bn A P from DE with respect to the feasible region; in which the parameters are linearly independent (Asprey and Macchietto,

A. Holmqvist et al. / Chemical Engineering Science 94 (2013) 316–329

2000; Vajda et al., 1989): n

Nb

P ¼ fb A R

Nb

: rank Jðb Þ ¼ Nb 4:(b A R , UðbÞ r Uðb Þg n

n

ð8Þ

the gradient-based algorithm lsqcurvefit in the Optimization TM for MATLABs (The MathWorks Inc, 2010b) was utiToolbox lized to yield the final least-squares estimates, b^ .

5. Results and discussion 5.1. Parameter estimation and statistical model validation Parameter estimation of the kinetic parameters (Eq. (1)) involved in the ALD gas–surface reactions was performed with Table 1 Outline of the experimental design. The ALD pulse sequence considered for the growth of ZnO was the following ðDt ZnðC2 H5 Þ2 , Dt N2 , Dt H2 O , Dt N2 Þ including a nitrogen purge

DtN2 ¼ 4:0  101 s between the non-overlapping alternate injections of the precursors for all cases E. All experiments were performed with the process operating parameters: /pS9Gout ¼ 2:0  102 Pa, ½Q_ ZnðC2 H5 Þ2 , Q_ H2 O , Q_ N2  ¼ ½0:11,0:16,5:0  102 sccm [standard cubic centimeters per minute at STP] (cf. Holmqvist et al., in press for physical interpretation). The uniformity factors for the experimentally observed response, UF y^ , and the model generated response, UFy, and the relative model error, dy, for all calibration sets are also presented.

DtZnðC2 H5 Þ2  101

DtH2 O  101 (s)

T  10  2 (1C)

UF y^ (%)

UFy (%)

dy

(s) 1 2.0 2a 4.0 3 6.0 4 12.0 5 18.0 6 6.0 7 6.0 8 6.0

4.0 4.0 4.0 4.0 4.0 1.0 2.0 4.0

1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0

97.1 94.7 93.8 90.1 94.7 92.3 93.7 93.8

94.4 95.6 95.9 96.4 95.9 96.2 95.8 95.6

5.9 4.1 6.7 5.7 2.7 9.9 9.3 6.1

9 1.0 10 2.0 11a 4.0 12 6.0 13 12.0 14 6.0 15 6.0 16 6.0

4.0 4.0 4.0 4.0 4.0 1.0 2.0 4.0

1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5

96.8 95.7 98.4 98.3 99.0 97.8 95.7 98.4

96.9 97.8 99.1 99.4 99.9 95.8 99.5 99.4

2.2 5.3 1.9 2.2 0.76 3.0 7.3 2.2

17 1.0 18 2.0 19 4.0 20 6.0 21 12.0 22 6.0 23a 6.0 24 6.0

4.0 4.0 4.0 4.0 4.0 1.0 2.0 4.0

2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0

97.7 97.5 96.8 98.9 99.1 99.0 99.1 99.0

97.9 98.4 99.1 99.0 99.6 99.0 99.6 99.1

2.8 2.9 5.9 3.8 5.1 1.0 1.3 3.8

E

a

(%)

Experimental validation set.

319

the hybrid numerical procedure described in Section 4. The experimental calibration data sets, y^ E , 8E 2 = f2; 11,23g used for evaluation of FðbÞ are reported in Holmqvist et al. (in press) and the operating conditions associated with the E th index are listed in Table 1. The remaining data sets are used for experimental model validation (see Section 5.1.3). 5.1.1. Assessment of the accuracy and reliability of the parameter estimates The statistical significance and reliability of b^ is assessed by evaluating the confidence statements presented in Section 3.1. A large confidence interval (Eq. (5)) or region (Eqs. (4) and (6)) suggests either that a parameter estimate cannot be discriminated from another with the experimental design considered or that there is a high uncertainty in the precision associated to this parameter, i.e. the model is inadequate in distinguishing a phenomenon in the system (Bard, 1974; Bates and Watts, 1988). The approximate marginal confidence intervals (Eq. (5)) and the corresponding margin of error, eb^ ¼ Sıı ðb^ Þ1=2 t ðj;1a=2Þ , with a sigı nificance level of a ¼ 5:0  102 and 8ıA N b are reported in Table 2. The narrow confidence intervals for parameters b^ ı and 8ıA f1; 2, . . . ,5g indicate that these estimates are adequately determined. Consistently larger relative errors, though still remaining moderate, were obtained for parameters b^ ı and 8ıA f6; 7, . . . ,9g, i.e. the kinetic parameters governing the rehydroxylation reaction (see e.g. Holmqvist et al., in press; Deminsky et al., 2004; Matero et al., 2000; Rahtu et al., 2001), and these are thus the most difficult parameters to statistically discriminate against others satisfactorily with the experimental design considered. In addition, the maximum molar concentration of surface sites L can be deduced from the physical properties determined from the analysis of film growth experiments (Coltrin et al., 2008). The experimental surface site density, determined from the density of deposited ZnO films, rs ¼ 5:4  103 kg m3 at T ¼423.2 K in ~ 2=3 N~ ¼ ¨ Torndahl et al. (2007), was found to be L ¼ ðrs M 1 s NÞ 5 2 1:94  10 mol m . This value shows reasonable agreement with the final estimate b^ 1 presented in Table 2, further verifying that this estimate was adequately determined. Another aspect affecting the reliability of b^ is parameter correlation, which may originate from an inappropriate model representation and/or experimental design (Schwaab and Pinto, 2007). Despite the implication of the intrinsic mathematical structure of the Arrhenius equation (Nagy and Tura´nyi, 2011; Varga et al., 2011), the majority of the elements of the parameter correlation matrix, C ı‘ ðb^ Þ and 8ıa ‘ A Nb , are moderately low (see P Table 2) and yielding an overall correlation level C ¼ ð ı a ‘ C ı‘ ðb^ Þ2 ðN2b N b Þ1 Þ1=2 (Pritchard and Bacon, 1978) of 0.43. This implies that the characteristically high parameter correlations, almost equal to one (He´berger et al., 1987), were successfully

Table 2 Regression analysis of the parameter least-squares estimates, b^ . The correlation matrix Cðb^ Þ and normalized margin of error e b^ were determined at a significance level of

a ¼ 5  102 . The specific reaction rate kT ref,i ,i and 8i A f1; 2, . . . ,4g are defined at the reference temperature T ref,i ¼ 4:19  102 K. According to Elam and George (2003), the variable numbers of surface hydroxyl groups that react with the ZnðC2 H5 Þ2 adsorptive precursor molecule in the elementary surface reactions (cf. Holmqvist et al., in press) was defined as n ¼ 1:37. Parameter units: kT ref,1 ,1 ððmol m2 Þ1n Pa1 s1 Þ; kT ref,2 ,2 ððmol m2 Þn1 Pa1 s1 Þ; kT ref,3 ,3 ððmol m2 Þ1 s1 Þ; kT ref,4 ,4 ðPa1 s1 Þ. ı

Parameter b

Parameter estimates b^

Margin of error e b^ (%)

Approximate correlation matrix Cðb^ Þ

1 2 3 4 5 6 7 8 9

L

1.53  10  5 8.16  10  6 3.57  104 1.33  10  6 4.04  104 3.45  10  6 2.67  104 4.68  10  7 3.79  104

0.183 1.24 0.172 1.16 0.155 3.29 2.56 5.21 4.95

1.00

kref,1 E1 kref,2 E2 kref,3 E3 kref,4 E4

0.63 1.00

0.74  0.53 1.00

0.46 0.12 0.21 1.00

0.52  0.22  0.15 0.52 1.00

 0.80 0.19  0.72  0.28  0.11 1.00

0.61 0.11 0.19 0.13 0.09  0.56 1.00

0.39 0.10 0.21  0.97 0.20  0.31 0.38 1.00

 0.35  0.06  0.23  0.30  0.17 0.36 0.29  0.59 1.00

320

A. Holmqvist et al. / Chemical Engineering Science 94 (2013) 316–329

decreased through reparameterization of the Arrhenius equation and definition of a reference temperature, see Holmqvist et al. (in press), Schwaab et al. (2008b), and Schwaab and Pinto (2007, 2008), leading to significant improvement in the precision of individual parameter estimates. In this study, T ref,i , 8i, were assigned the inverse average temperature over the experimental range analyzed (cf. Table 1). In order to understand the behavior of the relative errors and the parameter correlations, a more rigorous statistical evaluation of the parameter estimates is obtained by considering confidence regions (Eqs. (4) and (6)), described by function evaluations performed by the heuristic optimization method along the search, as discussed in Schwaab et al. (2008a). The likelihood 100ð1aÞ% joint confidence regions, with a significance level of a ¼ 5  102 , for pairs of parameters bı and 8ıA f1; 2,3; 6g are depicted in Fig. 1. From the figure it is evident that the significant reduction in parameter correlation, through the reparameterization method applied in Holmqvist et al. (in press), also provides an improvement in the elliptical representation of the confidence regions (Bates and Watts, 1988) associated with the pair of parameters b2 b3 (of the gas–surface reaction i¼1 Holmqvist et al., in press). Furthermore, it is noteworthy that the pairwise plots of the likelihood joint confidence regions clearly do not extend into negative parameter regions, which is required for interpreting physically meaningful parameters (Beale, 1960). However, when considering the 95% likelihood joint confidence region (Eq. (6)) associated with ıA f4; 8g depicted in Fig. 2, it is evident that the approximate 95% joint confidence region

Fig. 2. Nominal 95% likelihood confidence region for the normalized and centered parameters b4 and b8 . (– –) Approximate marginal 95% confidence interval ð7 e b^ Þ; (–) approximate joint 95% confidence region; (þ ) Normalized and ı centered least-squares estimates, b^ ı .

Fig. 1. Nominal 95% likelihood confidence region for pairs of normalized and centered parameters bı and 8ıA f1; 2,3; 6g; (– –) approximate marginal 95% confidence interval ð 7 e b^ Þ; (þ ) normalized and centered least-squares estimates, b^ ı . ı

A. Holmqvist et al. / Chemical Engineering Science 94 (2013) 316–329

(Eq. (4)) provides a poor representation, and thus, lead to erroneous conclusions about the confidence on model parameters. In Fig. 2, the point estimate of model parameters is placed at the center of the hyper-ellipsoid and the likelihood confidence region is non-convex and unbounded, since the parameter b8 does not have an upper limit. Hence, the model is unable to discriminate between high values of b8 , as high values of this parameter lead to a constant model response, y. The occurrence of the elongated likelihood joint confidence region reflects the parameter correlation (Emig and Hosten, 1974). In order to further reduce the parameter correlation and simultaneously improve the precision of the parameter estimates, two essential points are: (i) use an experimental design with extended information that enables discrimination (Franceschini and Macchietto, 2008) of the aforementioned parameters and (ii) perform the numerical optimization procedure proposed in Schwaab et al. (2008b) and Schwaab and Pinto (2007, 2008) to provide the optimum reference temperatures, T ref,i , for every ith Arrhenius reaction rate equation, respectively. An optimal experimental design and a set of optimal reference values holds the potential to yield precise estimates of uncorrelated parameters, and the likelihood joint confidence region would tend towards being elliptical. Moreover, on the subject of improving information content, in ^ to a statistical sense, obtained from the measured response, y, maximize the calibration parameter discrimination capacity, the ITRS guidelines promote the use of in situ diagnostics, including quartz crystal microbalance (QCM) (Aarik et al., 2006; Elam and George, 2003) and quadrupole mass spectroscopy (QMS) (Lei et al., 2006; Rahtu and Ritala, 2002). These monitoring techniques provide extensive dynamic information during the ALD pulse sequence, thereby improving the information content and enhancing the reliability of the parameter estimates. However, the main objective of the derived ALD process model in this article series was to accurately predict the response of the substrate spatially dependent film thickness profile to the variations of all the measured equipment-scale process operating parameters. For this reason, the ex situ instrumentation used for film characterization in Holmqvist et al. (in press) are superior to the aforementioned in situ diagnostics, since the latter measurement response provides no information regarding the nature of the spatial dependence of the deposition rate. Consequently, to further enhance the reliability of the parameter estimates, b^ , while still maintaining the precision of the spatial model response, the utilization of both in situ and ex situ diagnostic methods may become necessary.

5.1.2. Assessment of the accuracy and reliability of the predicted model response In order to verify the precision of the model response, the statistical uncertainties and their effects on the model response were addressed by considering the model 100ð1aÞ% approximate confidence band (Eq. (7)). The model generated spatially dependent film thickness profile, hs ðzÞ, for the least-squares estimates (see Table 2) and the associated 95% approximate confidence band are depicted in Fig. 3 along with the experimentally observed film ^ f^ , sampled at three spatially equidistant positions, thickness, y9 f^ ¼ ½16 , 36 , 56  9Gsub 9 (m), on the substrate, Gsub . Essentially, to statistically assess the reliability of the spatially dependent film thickness profile hs ðzÞ, it is necessary that the confidence bandwidth is moderately narrow and that it replicates the spatial dependency of the model response. According to Fig. 3 it is possible to statistically assure that the model output predicts a non-uniform film thickness profile for E ¼ 17, in contrast to the case of E ¼ 10. Hence, the accuracy and predictability of the model can thus be determined by analyzing the bandwidth and shape of the approximate inference band. In addition, the effect of the

321

uncertainty on the precision of the parameter estimates given by the model can be evaluated. The graphic representation of some of the experimental calibration data sets E A f1; 3,8; 9,10; 14,17; 20,24g in Fig. 3a–c indicates that the model output and the experimental data show good agreement. The comparison between the uniformity factors for the experimentally observed, UF y^ , and model generated response, UF y , in Table 1 also shows satisfactory agreement for the majority of the calibration data sets. The worst fit is found for the calibration data sets employing a temperature of 1:0  102 1C and short pulse duration of the adsorptive precursors, which give film profiles with the lowest uniformity factors. This is evident when visualizing the model accuracy, the associated bandwidth of the inference bands and when analyzing the relative error dyE and 8EA NE for each calibration index (see Table 1). On the other hand, the calibration sets with longer adsorptive precursor pulse duration, especially for the H2O precursor, have significantly lower relative errors, demonstrating that the model describes this part of the calibration region better. When growing oxide films by ALD in a viscous flow reactor with alternating injection of adsorptive precursors, Groner et al. (2004) reported that complete purging of reactants, especially H2O, takes significant longer time at low temperatures ð r1:0  102 1CÞ because of slower desorption from the reactor walls. Hence, purge times must be increased with decreasing temperature in order to prevent simultaneous presence of both precursors in the gas phase which results in regular CVD growth, comprising a higher growth rate than in true ALD (Elers et al., 2006; Jur and Parsons, 2011). In this study, the N2 purge time was held constant at Dt N2 ¼ 4:0  101 s (see Table 1) in the entire temperature region and local CVD growth may occur on the leading edge of the substrate in the cross-flow reactor, causing the sharp change in growth rate over the substrate, and thus degrading the overall uniformity. Consequently, this phenomena could explain the discrepancy between the model predicted response and the experimental measurement data at T ¼ 1:0  102 1C. A more comprehensive presentation of the reliability and the adequacy of the model response can be achieved by examining the distribution of the residual errors for the calibrated model. This can be graphically visualized by plotting the relative error dyE and 8EA N E for all calibration sets against the cumulative distribution, determined using the empirical cumulative distribution TM function ecdf in the Statistics Toolbox for MATLABs (The MathWorks Inc, 2010c). On the basis of the cumulative distribution of residual errors in Fig. 4, it is evident that 90% of the residual errors is below 7.3% and the only errors above this arise from experimental indices E A f6; 7g. However, the relative errors are below 10% for all experimental indices, confirming rather high predictability even for the worst cases studied. The fitted normal cumulative distribution function, using TM normcdf in the Statistics Toolbox for MATLABs (The MathWorks Inc, 2010c), depicted in Fig. 4, demonstrates reasonable agreement with the cumulative frequency of the relative model error, implying that the residual errors follow a normal distribution with an expected value and a variance (Hines et al., 2003), in accordance with the chi-square goodness-of-fit test as defined in Schwaab et al. (2006), Pinto et al. (2011), and Alberton et al. (2011). Thus, the underlying statistical assumptions on which the use of the least squares was based are not violated.

5.1.3. Experimental model validation The ultimate test of the adequacy of the model is to compare it to the validation set, and if the experimentally observed and model predicted response show good agreement, the model can be considered valid. The region of validity can subsequently be

322

A. Holmqvist et al. / Chemical Engineering Science 94 (2013) 316–329

1.14

1.49

1.48

1.03

1.34

1.34

0.93

1.2

1.2

0.82

1.06

1.07

0.71

0.92

0.93

1.98

2.27

2.22

1.86

2.1

2.06

1.74

1.93

1.9

1.62

1.76

1.74

1.49

1.59

1.58

2.09

2.21

2.21

1.96

2.07

2.07

1.84

1.92

1.92

1.72

1.78

1.78

1.59

1.63

1.64

1.26

2.27

2.01

1.16

2.17

1.93

1.06

2.07

1.86

0.95

1.97

1.79

0.85

1.87

1.72

Fig. 3. Spatially dependent film thickness profile hs ðzÞ as a function of the substrate boundary Gsub local coordinate variable z A ½0,9Gsub 9 (m). The operating conditions associated with the E th index are specified in Table 1. Model response hs ðzÞ (A˚ cycle  1) for: (a) calibration data sets performed at T ¼ 1:0  102 1C; (b) calibration data sets performed at T ¼ 1:5  102 1C; (c) calibration data sets performed at T ¼ 2:0  102 1C; (d) validation data sets. (þ ) Experimental response hs ðf^ Þ (A˚ cycle  1) at f^ ¼ ½16 , 36 , 56  9Gsub 9 (m); (–) spatial film thickness profile; (– –) 95% approximate confidence band (Eq. (7)); (––) uncertainty band for the expected response determined by uniformity distribution with the parameter bounds specified according to the marginal confidence intervals given in Table 2.

defined as the union between of the calibration region and the region covered by the validation experiments (Brereton, 2003). Fig. 3d visualizes the model response for the experimental validation data set with EA f2; 11,23g (see Table 1). The uncertainties enclosed bands were constructed from the extreme model response generated from Latin hypercube sampling (LHS) (McKay TM et al., 1979), using lhsdesign in the Statistics Toolbox for s MATLAB (The MathWorks Inc, 2010c), with the parameter boundaries specified in accordance with the marginal confidence intervals (see Eq. (5) and Table 2). The graphical presentation of the model response in Fig. 3d demonstrates that the performance of the model and the resulting predictability are satisfactory, with low relative residual errors, dy, and conformal uniformity factors for the model generated, UF y , and the experimentally observed, UF y^ , response, which confirms that the model is valid within the calibration region.

5.2. Mechanistic analysis of the ALD process model This section demonstrates the universal applicability of the mechanistic model of the continuous cross-flow ALD reactor (Suntola, 1992), derived in the first article of this series (Holmqvist et al., in press), with statistical reliable least-squares estimates (see Eq. (1) and Table 2) for the kinetic parameters involved in the elementary ALD gas–surface reactions. Based on the experimentally validated mathematical model, with adequate predictability of the spatial resolution of the film thickness profile (see Section 5.1), an extensive theoretical investigation was performed with the focus on mechanistic understanding of ALD process dynamics and the influence of equipment-scale process operating parameters on ALD film growth. The mechanistic model-based analysis complements available experimental results deduced from the literature, see e.g. Aarik et al. (2006),

A. Holmqvist et al. / Chemical Engineering Science 94 (2013) 316–329

2

100

80

60

40

20

0

2

4

6

8

10

12

1

1.8

0.9

1.6

0.8

1.4

0.7

1.2

0.6

1

0.5

0.8

0.4

0.6

0.3

0.4

0.2

0.2

0.1

0

0

323

0

50

100

150

200

250

0 300

14

Fig. 4. The cumulative distribution of the relative residual errors, dyE , and 8E A N E . (–) The cumulative frequency for the calibration data set; (– –) fitted normal TM cumulative distribution function normcdf in the Statistics Toolbox for MATLABs (The MathWorks Inc, 2010c).

Fig. 5. Dependence of the film mass increment on the process temperature during one ALD sequence. (–) The integral mean value of the mass increment /ms S over Gsub per ALD sequence and unit area; ð,Þ analogous to (–) but using the generalization of Fick’s law (cf. Holmqvist et al., in press); (– –) the degree of hydroxylation /gS (a.u.) after the H2 O exposure in the pulse sequence. Process operating parameters: /pS9Gout ¼ 3:0  102 Pa, /pS9Gin ¼ 4:3  102 Pa, ½Q_ ZnðC2 H5 Þ2 , Q_ H2 O , Q_ N2  ¼ ½0:10,0:15,5:0  102 sccm (standard cubic centimeters per minute at STP), pulse sequence, ½Dt ZnðC2 H5 Þ2 , Dt N2 , Dt H2 O , Dt N2  ¼ ½2:0,4:0,2:0,4:0  101 s.

Elam and George (2003), Jur and Parsons (2011), Kuse et al. (2003), Lei et al. (2006), and Rahtu et al. (2001).

5.2.1. General characteristics of the ALD film growth process In the context of deposition of ZnO films in continuous crossflow ALD reactors, the temperature, which is the most readily attribute to activation energies of the elementary reactions, as well as the process time parameters, i.e. precursor exposure and purge times, dependence have been undertaken seriously in Elam ¨ et al. (2002), Elam and George (2003), Torndahl et al. (2007), and Yousfi et al. (2000). According to Deminsky et al. (2004) and Matero et al. (2000), either the dehydroxylation reaction or the rehydroxylation reaction (surface reaction index i A f3; 4g, Holmqvist et al., in press) determines the overall growth rate in the high-temperature region. This is evident when considering the temperature dependence of the degree of hydroxylation, /gS, after the H2O sequence t A ½0, Dt ZnðC2 H5 Þ2 þ Dt N2 þ Dt H2 O , graphically presented in Fig. 5, and expressed as Z Z 1 yOH /gS ¼ @G@t ð9Þ 9Gsub 9 t Gsub yOH þ yO At the highest temperatures, dehydroxylation is so extensive that it overrides the increase in hydroxyl groups during the water sequence. In contrast, the overall growth rate decreases in the low-temperature region, which is attributed to the activation energy of the elementary surface reactions i A f1; 2g. A graphic representation of the transient production and consumption of fractional surface species, yk and k A fOH, ZnðC2 H5 Þ2n ,Og, during one complete ALD pulse sequence ðDt ¼ Dt ZnðC2 H5 Þ2 þ 2Dt N2 þ Dt H2 O Þ sampled at three different temperatures is presented in Fig. 6. It is evident from this figure that the fractional surface coverage of elementary oxygen, yO , is the source of the decrease in yOH at elevated temperatures. Because of the sequential nature of the surface reactions, the effect of the net yO production is twofold: the fractional surface coverage of yOH is reduced, which limits the number of adsorption sites for the adsorptive ZnðC2 H5 Þ2 precursor, thereby preventing the production of yZnðC2 H5 Þ2n , which successively reduces the number of

adsorption sites for the adsorptive H2O precursor in the subsequent pulse sequence. The pseudo-steady-state process corresponding to one pulse sequence can be seen in the figure, where the initial surface fraction composition is also governed at the end of the same pulse sequence. This is a consequence of the initial normally hydroxylated surface, resulting in a constant mass gain per cycle (MGPC) (Burton et al., 2009), and consistent with the ex situ XRR analysis in Holmqvist et al. (in press). However, all states governing the system’s behavior do not remain constant throughout the duration of the pulse sequence; therefore, the system is evidently not time-invariant. The response of the transient system is apparent when analyzing the mean mass deposited per unit area resolved at the level of each cycle, see Fig. 7. The theoretical trajectory resembles previous experimental QCM measurements for ZnO: i.e. a mass increment due to irreversible adsorption of the adsorptive Zn(C2H5)2 precursor, a flat portion of the trajectory during the purge exposure, and a decrease in mass during the adsorptive H2O precursor sequence (Elam and George, 2003; Yousfi et al., 2000). The general formalism of the sequential elementary surface reactions i A f1; 2g (Holmqvist et al., in press) incorporating a variable number of reacting hydroxyl groups, n, enables the surface chemistry to be analyzed during one pulse sequence, and the discrepancy between different reaction mechanisms to be investigated on the basis of mass increment during the respective half-reactions. The theoretical mass increment trajectories shown in Fig. 7 demonstrate the characteristic dependence of n A 0, nL ½, where nL represents the total number of ligands of the adsorptive organometallic precursor, and corresponding DM 2 DM 1 ratios 1 (Elam and George, 2003). Thus, the inherent structure of the model supports the utilization of the entire information content of observed mass increment trajectories when extracting detailed surface chemistry data, i.e. when determining n and the calibration parameter vector, b (Eq. (1)). An additional method of assessing n considers temporal and spatial integration over an arbitrary boundary G, of the reaction

324

A. Holmqvist et al. / Chemical Engineering Science 94 (2013) 316–329

100 1 0.8

80

0.6 0.4

60

0.2 0

40

1

20

0.8 0.6

0

0.4

15

0.2

12

0

9 6

1

3

0.8

0

0.6 0.4

0

0.2

0.4

0.8

1

1.2

0.2 0 0

0.2

0.6

0.8

1.2

Fig. 6. Integral mean value of the fractional surface coverage /yk S (a.u.) at Gsub during one ALD sequence at: (a) T ¼ 1:5  102 1C; (b) T ¼ 2:0  102 1C; (c) T ¼ 2:5  102 1C. (-) /yOH S; ð,Þ analogous to (-) but using the generalization of Fick’s law (cf. Holmqvist et al., in press); (- -) /yZnðC2 H5 Þ2n S; ð  Þ /yO S. The dotted lines perpendicular to the t-axis represent the individual adsorptive precursors and purge exposures of one pulse sequence. The process operating parameters are analogous to those defined in Fig. 5.

byproduct released during the surface reaction i¼1 (Holmqvist et al., in press) and the entire pulse sequence. Consequently, this quantity is governed by

n ¼ nL

Y Z Z W¼1

tW

G

pC2 H6 @G@t

ð1ÞW þ 1 ð10Þ

where W ¼ f1; 2g and t 1 A ½0, Dt ZnðC2 H5 Þ2 þ Dt N2  for the nominator respective t 2 A ½0, Dt ZnðC2 H5 Þ2 þ 2Dt N2 þ Dt H2 O  denominator temporal integrals. The trajectories of the integral mean value of the gas-phase composition at the outlet boundary, Gout , during one pulse sequence are presented graphically in Fig. 8. It is clear that the /pC2 H6 S9Gout trajectory is explicitly governed by xa,i /r i S9G sub and Eq. (10) provides the specified nominal value of n. The aforementioned method, which considers spatial integration, provides quantitative, chemically specific information and can be used as a supplementary theoretical analysis quantity in in situ QMS diagnostics (Lei et al., 2006; Rahtu et al., 2001). It is noteworthy that the spatial boundary integration or point evaluation can be utilized anywhere downstream of the substrate to adequately reproduce the output of QMS sampling systems.

Fig. 7. (a) Transient analysis of the integral mean mass increment /ms Sðng cm2 Þ at Gsub during one ALD sequence. (–) /ms S with the number of reacting hydroxyl groups n ¼ 1:37; (– –) /ms S with the number of reacting hydroxyl groups n ¼ 1:0. (b) Transient analysis of the integral mean gas-phase composition, /pa S (Pa), of the adsorptive precursor at Gin . (–) /pZnðC2 H5 Þ2 S; ð,Þ analogous to (–) but using the generalization of Fick’s law (cf. Holmqvist et al., in press); (– –) /pH2 O S. The dotted lines perpendicular to the t-axis represent the individual adsorptive precursors and purge exposures of one pulse sequence. The process operating parameters are analogous to those defined in Fig. 5 and T ¼ 1:5  102 1C.

Together, Figs. 7 and 8 reveal rigorous process dynamics of the ZnO ALD process that are essential process analysis, design and optimization. The comparison of the model response in Figs. 5–8, using the binary and multicomponent formulations (cf. Holmqvist et al., in press), respectively, for the mass flux show only minor deviations. This is due to the low molar composition of C2 H6 relative the remaining species, and hence, the rigorous, multicomponent formulation may be omitted, to a very good approximation, and one-way solute–solvent interactions suffice. Therefore, the discussion regarding the sensitivity of the process operating parameters in Section 5.2.2 will be addressed solely on the basis of the binary formulation for the mass flux.

5.2.2. Sensitivity of process operating parameters on ALD film growth In their experimental investigation of deposition efficiency in a continuous flow ALD reactor, Aarik et al. (2006) and Jur and Parsons (2011) observed the influence of carrier gas flow parameters on ALD growth rate. To verify whether the predicted model response resemble the experimental observations, a comprehensive univariate sensitivity analysis of the measured process operating parameters u ¼ ½/pS9Gout ,T, Q_ a , Dt a y , a A fN2 ,ZnðC2 H5 Þ2 ,H2 Og,

A. Holmqvist et al. / Chemical Engineering Science 94 (2013) 316–329

10

8

6

4

2

0 10 8 6 4 2 0

0

0.2

0.4

0.8

1

1.2

Fig. 8. (a) Transient analysis of the integral mean gas-phase composition, /pa S (Pa), at Gout during one ALD sequence. (–) /pZnðC2 H5 Þ2 S; ð,Þ analogous to (–) but using the generalization of Fick’s law (cf. Holmqvist et al., in press); (––) /pC2 H6 S  5; (–) /pH2 O S. (b) Transient analysis of the integral mean value of elementary surface reactions /r i S ðmol m2 s1 Þ at Gsub . (–) /r 1 S; (– –) /r2 S. The dotted lines perpendicular to the t-axis represent the individual adsorptive precursors and purge exposures of one pulse sequence. The process operating parameters are analogous to those defined in Fig. 5 and T ¼ 1:5  102 1C.

on the increment of film mass per unit area, resolved at the level of a single pulse sequence, /ms S was performed. Furthermore, the model-based analysis is extended to comprehend the mechanistic parameter impact on the spatially dependent film thickness profile, conveniently represented by the uniformity factor, UF. In Aarik et al. (2006), the influence of carrier gas mass flow, Q_ N2 , was evaluated on the basis of the linear flow rate, and in this study, expressed mathematically as the integral mean value of the velocity field, /vS9O , in the computational domain O (cf. Holmqvist et al., in press): Z

Z

1 X /vS9O ¼ v@O ¼ r Q_ a @O 9O9 O 9O9 O rAO a ¼ 1 STP, a 1

1

ð11Þ

where AO is the cross section area of O and r ¼ pMðRTÞ1 the local density. The predicted model response depicted in Fig. 9a shows that /ms S decreases, at an applied carrier gas pressure, /pS9Gout , at Gout , when increasing the linear flow rate of the carrier gas. Analogous dependencies were observed when increasing the carrier gas pressure while maintaining a constant linear flow rate, although in this case, the sensitivity of the carrier gas pressure on /ms S was weaker. This conclusion is consistent with the experimentally observed parametric behavior in Aarik et al. (2006).

325

However, it is reasonable to exclude the decrease in reactor 1 chamber residence time t ¼ ð9Gsub 9 þ29Gwall 9Þ/vS9O (s) at elevated linear flow rates as the underlying cause of the observed significant decrease in /ms S, since t 5 Dt a , a A fZnðC2 H5 Þ2 ,H2 Og, for the entire range of /vS9O studied. In contrast, the decrease in /ms S can be mechanistically correlated to the hydrodynamic interdependence governed by Eq. (11), which originates from the ideal gas law. At constant applied /pS9Gout and T, Eq. (11) states that the linear flow rate is directly proportional to the mass flow of the carrier gas, /vS9O pQ_ N2 . Consequently, the increase in Q_ N2 with /vS9O imposes a decrease in local mass fraction of the a th adsorptive precursors, oa , and ultimately resulting in reduced driving force in the gas–surface reaction rates, hence r i pðpca Þ9G and 8i a 3 (Holmqvist et al., in press). sub According to Eq. (11), a decrease in carrier gas pressure at Gout , at constant applied /vS9O and T, imposes a decrease in Q_ N2 and an increase in the overall local oa , respectively. Additionally, overall lower local p also implies that the binary gas diffusion coefficient is increased, since Dab pp1 (Bird et al., 1960), which increases the rates of impingement mass flux to the substrate surface, ja ¼ rDab roa (Holmqvist et al., in press). In summary, these mechanistic characteristics will promote a higher molar composition of the a th adsorptive precursors, ca , in the vicinity of the substrate, enhancing the rate of elementary surface reactions. In contrast, a decrease in the local carrier gas pressure will limit this driving force. However, the net increase in mass deposited, /ms S, can be mechanistically correlated to the inverse pressure dependence of the impingement mass flux, ja , which compensates for the reduction in reaction driving force, r i pðpca Þ9G for 8i a 3. sub The temperature dependence on /ms S is thoroughly discussed in Section 5.2.1. The interdependence of the process temperature on the applied pressure and linear flow rate of the carrier gas is graphically presented in Fig. 9b and c. It is evident that the sensitivity to /vS9O is higher than to /pS9Gout , and the absolute effect on /ms S is more pronounced at temperatures corresponding to high growth rates. Provided that /pS9Gout , /vS9O , and T are held constant, there are two main ways of varying the half-cycle average exposure dose for the a th adsorptive precursor (Adomaitis, 2010; Gordon et al., 2003; Jur and Parsons, 2011): Z Z 1 p @G@t ð12Þ /da S ¼ 9Gsub 9 Dta Gsub a either by changing the mass flow of the precursors, Q_ a , thereby changing the partial pressure, or by changing the duration of the pulse, Dt a . Hence, under prevailing operating conditions, the required minimum substrate saturation exposure dose of the a th adsorptive precursor is governed by Z Z 1 n  ðroa v þ ja Þ L¼ @G@t ð13Þ Ma 9Gsub 9 Dtsat, a Gsub As emphasized by Kuse et al. (2003), a net effect of Q_ a and Dt a on /ms S can only be observed provided that the saturation exposure dose according to Eq. (13) is not achieved. The interdependence of mass flow and the pulse duration of ZnðC2 H5 Þ2 adsorptive precursor on /ms S is graphically presented in Fig. 9d. It is evident that /ms S reached a plateau as saturation of the surface sites is achieved, which demonstrates the characteristic self-limiting nature of the ALD processes (Lim et al., 2003; Puurunen, 2005). As emphasized by Cleveland et al. (2012), the across substrate uniformity is a useful process metric to evaluate the self-limiting behavior of the ALD process. In agreement with experimental observations in Henn-Lecordier et al. (2011), enhanced uniformity at higher precursor exposure doses, /dZnðC2 H5 Þ2 S, were obtained in Fig. 10d. However, it is noteworthy that the film thickness

326

A. Holmqvist et al. / Chemical Engineering Science 94 (2013) 316–329

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0 50

0 50

180

310

440

570

700

70

59

48

37

26

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0 15

0 1.5

26

37

48

59

70

250

210

170

130

90

180

15

50

1.2

310

0.9

440

0.6

570

0.3

700

0

250

0

90

210

130

50

170

16

4

12

20

8

Fig. 9. Normalized integral mean mass increment /ms S (a.u.) as a function of (a) integral mean value of the velocity field /vS9O ðm s1 Þ on O and pressure at the outlet /pS9Gout (Pa); (b) process temperature T (1C) and outlet pressure /pS9Gout (Pa); (c) process temperature T (1C) and integral mean value of the velocity field /vS9O ðm s1 Þ on O; (d) inlet mass flow Q_ ZnðC2 H5 Þ2 (sccm) and pulse duration DtZnðC2 H5 Þ2 (s). The values of the parameters not varied in each case were assigned according to the reference conditions: temperature T ¼ 1:5  102 1C, pressure at the outlet /pS9Gout ¼ 3:0  102 Pa, integral mean value of the velocity field /vS9O ¼ 3:0  101 m s1 on O, inlet adsorptive precursor mass flow ½Q_ ZnðC2 H5 Þ2 , Q_ H2 O  ¼ ½5:0,5:0 sccm, pulse sequence ½Dt ZnðC2 H5 Þ2 , Dt N2 , Dt H2 O , Dt N2  ¼ ½2:0,4:0,2:0,4:0  101 s).

uniformity factor, UF, is not strictly increasing with the precursor exposure dose. Furthermore, the interdependence of /pS9Gout and /vS9O depicted in Fig. 10a reveals that UF approaches unity at the lower parameter boundaries. By analysing the spatially dependent film thickness profiles depicted in Fig. 3, it is evident that thickness uniformity enhancement imposes higher reaction driving forces, r i pðpca Þ9G and 8i a 3, as z-9Gsub 9 (m). As aforementioned, the sub decrease in /pS9Gout and /vS9O promoted enhanced mass flux, ja , to the substrate surface and mass fraction, oa , of the a th precursor, respectively, and thus, explains the parametric dependence on UF. In addition, the UF temperature dependence resembles that of /ms S (cf. Figs. 9 and 10), however the impact is not that pronounced in the lower temperature region, in which the overall growth rate is attributed to the activation energy of r i and i A f1; 2g (cf. Section 5.2.1).

6. Concluding remarks This paper demonstrates the universal applicability of the mechanistic model of the continuous cross-flow reactor system F-120 by ASM Microchemistry Ltd. Suntola (1992) derived in the first article of this series (Holmqvist et al., in press). According to the ITRS initiative, models of ALD processes must have a high level of accuracy if they are to be used to assist in process

analyses, design or optimization. For this reason, a systematic methodology to support the development and statistical validation of ALD process models was proposed. Within this framework, the statistical uncertainty of least-squares estimates of kinetic parameters involved in the ALD gas–surface reaction mechanism governing ZnO film growth, and its ultimate impact on the model response have been assessed. The comprehensive model assessment reported here showed that the expected high parameter correlation, originating from the complex interdependence between elementary surface reactions involved in ALD and the intrinsic mathematical structure of the Arrhenius equation, was successfully diminished through reparameterization of the Arrhenius equation, and hence, leading to the improvement in precision of the parameter estimates. It is also shown that the model predictions of the spatially dependent film thickness profile are in good agreement with both calibration and validation experimental data, respectively, under a wide range of operating conditions. The novel model-based methodology presented in this paper constitutes a general platform for mechanistic investigations of the detailed process dynamics of a continuous cross-flow ALD reactor. The extensive theoretical analysis exhibits that the experimentally validated model, successfully reproduces the detailed process dynamics revealed by modern in situ monitoring technologies, as advocated by the ITRS initiative, and thereby provides a supplementary analysis tool. Furthermore, the univariate sensitivity analysis

A. Holmqvist et al. / Chemical Engineering Science 94 (2013) 316–329

a

b

1

1

0.96

0.96

0.92

0.92

0.88

0.88

0.84

0.84

0.8 50

0.8 50

180

310

440

570

700

15

26

37

48

59

180

70

c

327

310

440

570

700

50

210

90

170

250

130

4

16

8

0

12

d 1

1

0.96

0.96

0.92

0.92

0.88

0.88

0.84

0.84

0.8 15

0.8 1.5

26

37

48

59

70

50

90

130

170

210

1.2

250

0.9

0.6

0.3

0

20

Fig. 10. Across substrate film thickness uniformity, UF, as a function of (a) integral mean value of the velocity field /vS9O ðm s1 Þ on O and pressure at the outlet /pS9Gout (Pa); (b) process temperature T (1C) and outlet pressure /pS9Gout (Pa); (c) process temperature T (1C) and integral mean value of the velocity field /vS9O ðm s1 Þ on O; (d) inlet mass flow Q_ ZnðC2 H5 Þ2 (sccm) and pulse duration Dt ZnðC2 H5 Þ2 (s). The values of the parameters not varied in each case were assigned according to the reference condition in Fig. 9.

revealed the mechanistic dependence of the measured process operating parameters on ALD film growth performance metrics. Hence, the presented model-based framework also serves as a means to guide future research efforts in the field of ALD process optimization. On the basis of the characteristic timescales of the mass transport and the surface reaction mechanism considered, it is possible to determine the optimal process operating parameters to achieve optimal film growth rate and adsorptive precursor utilization, while simultaneously controlling film thickness uniformity across the substrate.

g hs J ja ki

response function film height parameter Jacobian matrix diffusive mass flux vector reaction rate constant

Ma

molar mass

ms N~

film mass increment Avogadro’s number

p Q_

Ai C Dab Ei F

frequency factor in the Arrhenius equation Nb  N b parameter correlation matrix binary diffusivity activation energy

system of differential algebraic equations F ðNb , j;aÞ Fisher’s F-distribution

a

R

Roman letters ðmol m2 Þ1nads,i Pa1 s1 , ðmol m2 Þ1ndes,i s1 –

s2 T t t ðj;a=2Þ UF

m2 s1 J mol – –

1

ðmol m2 Þ1nads,i Pa1 s1 , ðmol m2 Þ1ndes,i s1

ni Nomenclature

– m – kg m2 s1

u v W w

1

kg mol kg m2 1

mol surface reaction adsorption – and desorption order pressure Pa volumetric flow rate at STP Nm3 s1 1 universal gas constant J mol K1 the variance estimate of s2 temperature time Student’s t-distribution spatial film profile uniformity factor design variables velocity vector weight matrix algebraic variables

– K s – – – m s1 – –

328

x y

A. Holmqvist et al. / Chemical Engineering Science 94 (2013) 316–329

state variables measured state variables

– –

Greek letters

a b

g G da dy

eb^ ı z yk L

t n

xi

r R F

j ca O

oa

significance level calibration parameter vector degree of hydroxylation portioned boundary half-cycle average precursor dose relative error residual normalized margin of error substrate boundary Gsub local coordinate variable fractional surface coverage of surface species maximum molar concentration of surface sites residence time numbers of surface OH groups reacting with each ZnðC2 H5 Þ2 molecule surface reaction stoichiometric coefficient density parameter covariance matrix weighted sum of squared residuals degrees of freedom molar fraction of gaseous species spatial domain mass fraction of gaseous species

– – – – Langmuir – – m – mol m2

s –

– kg m3 – – – – – –

Subscripts and superscripts 0

a, b i ı E

k ref STP s

initial value gaseous species index surface reaction index calibration parameter index calibration and validation set index surface species index state at the reference temperature standard temperature and pressure solid

Acknowledgement This work has been supported by the Swedish Research Council under Grant no. 2006-3738. References ¨ Aarik, J., Aidla, A., Kasikov, A., Mandar, H., Rammula, R., Sammelselg, V., 2006. Influence of carrier gas pressure and flow rate on atomic layer deposition of HfO2 and ZrO2 thin films. Appl. Surf. Sci. 252 (16), 5723–5734.

Aarik, J., Siimon, H., 1994. Characterization of adsorption in flow type atomic layer epitaxy reactor. Appl. Surf. Sci. 81 (3), 281–287. Adomaitis, R.A., 2010. Development of a multiscale model for an atomic layer deposition process. J. Cryst. Growth 312 (8), 1449–1452. Alberton, A.L., Schwaab, M., Biscaia Jr., E.C., Pinto, J.C., 2010. Sequential experimental design based on multiobjective optimization procedures. Chem. Eng. Sci. 65 (20), 5482–5494. Alberton, A.L., Schwaab, M., Loba~ o, M.W.N., Pinto, J.C., 2011. Experimental design for the joint model discrimination and precise parameter estimation through information measures. Chem. Eng. Sci. 66 (9), 1940–1952. Asprey, S.P., Macchietto, S., 2000. Statistical tools for optimal dynamic model building. Comput. Chem. Eng. 24 (2–7), 1261–1267. Bard, Y., 1974. Nonlinear Parameter Estimation. Academic Press, New York. Bates, D.M., Watts, D.G., 1988. Nonlinear Regression Analysis and Its Applications. John Wiley & Sons, Inc., New York. Beale, E.M.L., 1960. Confidence regions in non-linear estimation. J. R. Stat. Soc. B 22 (1), 41–88. Bird, R.B., Stewart, W.E., Lightfoot, E.N., 1960. Transport Phenomena, 2nd ed. John Wiley & Sons, Inc., New York. Brereton, R.G., 2003. Chemometrics, Data Analysis for the Laboratory and Chemical Plant, 1st ed. John Wiley & Sons, Inc., Chichester. Burton, B.B., Goldstein, D.N., George, S.M., 2009. Atomic layer deposition of MgO using bis(ethylcyclopentadienyl)magnesium and H2O. J. Phys. Chem. C 113 (5), 1939–1946. Buzzi-Ferraris, G., Manenti, F., 2009. Kinetic models analysis. Chem. Eng. Sci. 64 (5), 1061–1074. Chappell, M.J., Gunn, R.N., 1998. A procedure for generating locally identifiable reparameterisations of unidentifiable non-linear systems by the similarity transformation approach. Math. Biosci. 148 (1), 21–41. Cleveland, E.R., Henn-Lecordier, L., Rubloff, G.W., 2012. Role of surface intermediates in enhanced, uniform growth rates of TiO2 atomic layer deposition thin films using titanium tetraisopropoxide and ozone. J. Vac. Sci. Technol. A 30 (1), 01A150. Coello, C.A.C., Lamont, G.B., Veldhuizen, D.A.V., 2007. Evolutionary Algorithms for Solving Multi-Objective Problems, 2nd ed. Springer, New York. Coltrin, M.E., Hsu, J.W.P., Scrymgeour, D.A., Creighton, J.R., Simmons, N.C., Matzke, C.M., 2008. Chemical kinetics and mass transport effects in solution-based selective-area growth of ZnO nanorods. J. Cryst. Growth 310 (3), 584–593. COMSOL AB, 2008. COMSOL Multiphysics User’s Guide, Version 3.5a. COMSOL AB, Tegne´rgatan 23, SE-111 40, Stockholm. Davis, M.E., 1984. Numerical Methods and Modeling for Chemical Engineers, 1st ed. John Wiley & Sons, Inc., New York. Deminsky, M., Knizhnik, A., Belov, I., Umanskii, S., Rykova, E., Bagatur’yants, A., et al., 2004. Mechanism and kinetics of thin zirconium and hafnium oxide film growth in an ALD reactor. Surf. Sci. 549 (1), 67–86. Draper, D., 1995. Assessment and propagation of model uncertainty (with discussion). J. R. Stat. Soc. B 57 (1), 45–97. Draper, N.R., Guttman, I., 1995. Confidence intervals versus regions. J. R. Stat. Soc. D 44 (3), 399–403. Elam, J.W., George, S.M., 2003. Growth of ZnO/Al2O3 alloy films using atomic layer deposition techniques. Chem. Mater. 15 (4), 1020–1028. Elam, J.W., Sechrist, Z.A., George, S.M., 2002. ZnO/Al2O3 nanolaminates fabricated by atomic layer deposition: growth and surface roughness measurements. Thin Solid Films 414 (1), 43–55. Elers, K.E., Blomberg, T., Peussa, M., Aitchison, B., Haukka, S., Marcus, S., 2006. Film uniformity in atomic layer deposition. Chem. Vapor Deposition 12 (1), 13–24. Elliott, S.D., 2007. Models for ald and mocvd growth of rare earth oxides. In: Fanciulli, M., Scarel, G. (Eds.), Rare Earth Oxide Thin Films. Springer, Berlin, Heidelberg, New York, pp. 73–86. Emig, G., Hosten, L.H., 1974. On the reliability of parameter estimates in a set of simultaneous nonlinear differential equations. Chem. Eng. Sci. 29 (2), 475–483. Franceschini, G., Macchietto, S., 2008. Model-based design of experiments for parameter precision: state of the art. Chem. Eng. Sci. 63 (19), 4846–4872. Gordon, R.G., Hausmann, D., Kim, E., Shepard, J., 2003. A kinetic model for step coverage by atomic layer deposition in narrow holes or trenches. Chem. Vapor Deposition 9 (2), 73–78. Groner, M.D., Fabreguette, F.H., Elam, J.W., George, S.M., 2004. Low-temperature Al2O3 atomic layer deposition. Chem. Mater. 16 (4), 639–645. Hamilton, D.C., Watts, D.G., Bates, D.M., 1982. Accounting for intrinsic nonlinearity in nonlinear regression parameter inference regions. Ann. Stat. 10 (2), 386–393. He´berger, K., Keme´ny, S., Vido´czy, T., 1987. On the errors of Arrhenius parameters and estimated rate constant values. Int. J. Chem. Kinet. 19 (3), 171–181. Henn-Lecordier, L., Anderle, M., Robertson, E., Rubloff, G.W., 2011. Impact of parasitic reactions on wafer-scale uniformity in water-based and ozone-based atomic layer deposition. J. Vac. Sci. Technol. A 29 (5), 051509. Hibbert, D.B., 1993. Genetic algorithms in chemistry. Chemometrics Intelligent Lab. Syst. 19 (3), 277–293. Hines, W.W., Montgomery, D.C., Goldsman, D.M., Borror, C.M., 2003. Probability and Statistics in Engineering, 4th ed. John Wiley & Sons, Inc., New York. ¨ ¨ Holmqvist, A., Torndahl, T., Stenstrom, S. A model-based methodology for the analysis and design of atomic layer deposition processes—part I: mechanistic modelling of continuous flow reactors. Chem. Eng. Sci, http://dx.doi.org/10. 1016/j.ces.2012.07.015, in press.

A. Holmqvist et al. / Chemical Engineering Science 94 (2013) 316–329

Jur, J.S., Parsons, G.N., 2011. Atomic layer deposition of Al2O3 and ZnO at atmospheric pressure in a flow tube reactor. ACS Appl. Mater. Interfaces 3 (2), 299–308. Kuse, R., Kundu, M., Yasuda, T., Miyata, N., Toriumi, A., 2003. Effect of precursor concentration in atomic layer deposition of Al2O3. J. Appl. Phys. 94 (10), 6411–6416. Lei, W., Henn-Lecordier, L., Anderle, M., Rubloff, G.W., Barozzi, M., Bersani, M., 2006. Real-time observation and optimization of tungsten atomic layer deposition process cycle. J. Vac. Sci. Technol. B 24 (2), 780–789. Lim, B.S., Rahtu, A., Gordon, R.G., 2003. Atomic layer deposition of transition metals. Nat. Mater. 2 (11), 749–754. Lim, J.-W., Park, J.-S., Kang, S.-W., 2000. Kinetic modeling of film growth rates of TiN films in atomic layer deposition. J. Appl. Phys. 87 (9), 4632–4634. ¨ M., Sajavaara, T., 2000. Effect of water dose Matero, R., Rahtu, A., Ritala, M., Leskela, on the atomic layer deposition rate of oxide thin films. Thin Solid Films 368 (1), 1–7. McKay, M.D., Beckman, R.J., Conover, W.J., 1979. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 21 (2), 239–245. Moles, C.G., Mendes, P., Banga, J.R., 2003. Parameter estimation in biochemical pathways: a comparison of global optimization methods. Genome Res. 13, 2467–2474. Nagy, T., Tura´nyi, T., 2011. Uncertainty of Arrhenius parameters. Int. J. Chem. Kinet. 43 (7), 359–378. Park, H.-S., Min, J.-S., Lim, J.-W., Kang, S.-W., 2000. Theoretical evaluation of film growth rate during atomic layer epitaxy. Appl. Surf. Sci. 158, 81–91. Pinto, J., Loba~ o, M., Alberton, A., Schwaab, M., Embiruu, M., Vieira De Melo, S., 2011. Critical analysis of kinetic modeling procedures. Int. J. Chem. React. Eng. 9 (A87), 1–33. Price, K.V., 1999. An introduction to differential evolution. In: Corne, D., Dorigo, M., Glover, F. (Eds.), New Ideas in Optimization. McGraw-Hill, London, pp. 79–108. Price, K.V., Storn, R.M., Lampinen, J.A., 2005. Differential Evolution: A Practical Approach to Global Optimization, 1st ed. Springer, Berlin. Pritchard, D.J., Bacon, D.W., 1978. Prospects for reducing correlations among parameter estimates in kinetic models. Chem. Eng. Sci. 33 (11), 1539–1543. Puurunen, R.L., 2005. Surface chemistry of atomic layer deposition: a case study for the trimethylaluminum/water process. J. Appl. Phys. 97 (12), 121301. Puurunen, R.L., Vandervorst, W., 2004. Island growth as a growth mode in atomic layer deposition: a phenomenological model. J. Appl. Phys. 96 (12), 7686–7695. Rahtu, A., Alaranta, T., Ritala, M., 2001. In situ quartz crystal microbalance and quadrupole mass spectrometry studies of atomic layer deposition of aluminum oxide from trimethylaluminum and water. Langmuir 17 (21), 6506–6509. Rahtu, A., Ritala, M., 2002. Reaction mechanism studies on the zirconium chloride– water atomic layer deposition process. J. Mater. Chem. 12 (5), 1484–1489. ¨ M., 2002. Handbook of Thin Film Materials. vol. 1. Academic Ritala, M., Leskela, Press, New York. Schiesser, W.E., 1991. The Numerical Method of Lines: Integration of Partial Differential Equations, 1st ed. Academic Press, San Diago. Schwaab, M., Biscaia, J.E.C., Monteiro, J.L., Pinto, J.C., 2008a. Nonlinear parameter estimation through particle swarm optimization. Chem. Eng. Sci. 63 (6), 1542–1552.

329

Schwaab, M., Lemos, L.P., Pinto, J.C., 2008b. Optimum reference temperature for reparameterization of the Arrhenius equation. Part 2: problems involving multiple reparameterizations. Chem. Eng. Sci. 63 (11), 2895–2906. Schwaab, M., Pinto, J.C., 2007. Optimum reference temperature for reparameterization of the Arrhenius equation. Part 1: problems involving one kinetic constant. Chem. Eng. Sci. 62 (10), 2750–2764. Schwaab, M., Pinto, J.C., 2008. Optimum reparameterization of power function models. Chem. Eng. Sci. 63 (18), 4631–4635. Schwaab, M., Silva, F.M., Queipo, C.A., Barreto Jr., A.G., Nele, M., Pinto, J.C., 2006. A new approach for sequential experimental design for model discrimination. Chem. Eng. Sci. 61 (17), 5791–5806. Semiconductor Industry Association, 2011. Modeling and simulations. In: The International Technology Roadmap for Semiconductors, 2011th ed. San Jose, CA, pp. 1–45, Avaliable online: /http://public.itrs.net/reports.htmlS. Siimon, H., Aarik, J., 1995. Modelling of precursor flow and deposition in atomic layer deposition reactor. J. Phys. IV 05 (C5), 245–252. Siimon, H., Aarik, J., 1997. Thickness profiles of thin films caused by secondary reactions in flow-type atomic layer deposition reactors. J. Phys. D 30 (12), 1725–1728. Storn, R., Price, K., 1997. Differential evolution—a simple and efficient heuristic for global optimization over continuous spaces. J. Global Optim. 11 (4), 341–359. Suntola, T., 1992. Atomic layer epitaxy. Thin Solid Films 216 (1), 84–89. The MathWorks, Inc, 2010a. MATLAB R2010b Documentation, Version 7.11. The MathWorks, Inc, 3 Apple Hill Drive, Natick, MA 01760-2098. TM The MathWorks, Inc, 2010b. Optimization Toolbox User’s Guide, Version 7.11. The MathWorks, Inc, 3 Apple Hill Drive, Natick, MA 01760-2098. TM The MathWorks, Inc, 2010c. Statistic Toolbox User’s Guide, Version 7.11. The MathWorks, Inc, 3 Apple Hill Drive, Natick, MA 01760-2098. ¨ ¨ Torndahl, T., Platzer-Bjorkman, C., Kessler, J., Edoff, M., 2007. Atomic layer deposition of Zn1  xMgxO buffer layers for Cu(In,Ga)Se2 solar cells. Prog. Photovoltaics: Res. Appl. 15 (3), 225–235. Vajda, S., Rabitz, H., Walter, E., Lecourtier, Y., 1989. Qualitative and quantitative identifiability analysis of nonlinear chemical kinetic-models. Chem. Eng. Commun. 83 (1), 191–219. Varga, L., Szabo´, B., Zse´ly, I., Zemple´ni, A., Tura´nyi, T., 2011. Numerical investigation of the uncertainty of Arrhenius parameters. J. Math. Chem. 49, 1798–1809. Watts, D.G., 1994. Estimating parameters in nonlinear rate equations. Can. J. Chem. Eng. 72 (4), 701–710. Witkowski, W.R., Allen, J.J., 1993. Approximation of parameter uncertainty in nonlinear optimization-based parameter estimation schemes. AIAA J. 31 (5), 947–950. Xu, K., Ye, P.D., 2010. Theoretical study of atomic layer deposition reaction mechanism and kinetics for aluminum oxide formation at graphene nanoribbon open edges. J. Phys. Chem. B 114 (23), 10505–10511. Ylilammi, M., 1995. Mass transport in atomic layer deposition carrier gas reactors. J. Electrochem. Soc. 142 (7), 2474–2479. Yousfi, E.B., Fouache, J., Lincot, D., 2000. Study of atomic layer epitaxy of zinc oxide by in-situ quartz crystal microgravimetry. Appl. Surf. Sci. 153 (4), 223–234. Zienkiewicz, O.C., Taylor, R.L., 2000a. The Finite Element Method—Fluid Dynamics, 5th ed., vol. 3. Butterworth-Heinemann, Oxford. Zienkiewicz, O.C., Taylor, R.L., 2000b. The Finite Element Method—The Basis, 5th ed., vol. 1. Butterworth-Heinemann, Oxford.