ABSTRACT
Title of dissertation: REACTION FACTORIZATION
FOR DISTRIBUTED DEPOSITION SYSTEMS:
APPLICATION TO COPPER FILM GROWTH
David Arana-Chavez, Doctor of Philosophy, 2015
Dissertation directed by: Professor Raymond Adomaitis
Department of Chemical and Biomolecular Engineering
A model reduction methodology based on a Gauss-Jordan reaction factoriza-
tion for thin-film deposition reaction systems is developed in this thesis. The fac-
torization generates a transformation matrix that is used to create a new coordinate
system that guides the separation of the deposition process time scales by decoupling
the net-forward reaction rates to the greatest extent possible. The new coordinate
space enables recasting the original model as a singular perturbation problem and
consequently as a semi-explicit system of differential-algebraic equations (DAE) for
the dominant dynamics in the pseudo-equilibrium limit. Additionally, the factor-
ization reveals conserved quantities in the new reaction coordinate system as well
as potential structural problems with the deposition reaction network.
The reaction factorization methodology is formulated to be suitable for appli-
cation to dynamic, spatially distributed reaction systems. The factorization provides
a rigorous pathway to decouple the time evolution and the spatial distributions of
deposition systems when the dynamics of reactor-scale gas-phase transport are fast
relative to the deposition process. Moreover, the factorization approach provides a
solution to the problem of formulating Danckwerts-type boundary conditions where
gas-phase equilibrium reactions are important.
The reaction factorization is used to study the chemical vapor deposition of
copper on a tubular hot-wall reactor using copper iodide as the Cu precursor. A
film-growth mechanism is proposed from experimental observations that the copper
films deposited on quartz substrates suggest a Volmer-Weber growth mode. A model
based on this mechanism is used to track spatial distribution of the average Cu island
size in the reactor.
The rate expressions used in the Cu deposition model are determined using
absolute rate theory. To carry out these calculations in an organized manner, a
library of object-oriented classes are created in the Python programming language.
REACTION FACTORIZATION FOR DISTRIBUTED
DEPOSITION SYSTEMS: APPLICATION TO COPPER FILM
GROWTH
by
David Arana-Chavez
Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
2015
Advisory Committee:
Professor Raymond A. Adomaitis, Chair/Advisor
Professor Panagiotis Dimitrakopoulos
Professor Jeffery Klauda
Professor Michael R. Zachariah
Professor Oded Rabin
c© Copyright by
David Arana-Chavez
2015
Acknowledgments
Above all I thank God for without Him nothing would be possible.
Infinite thanks to my advisor, Dr. Raymond Adomaitis (Ray), for the con-
tinuous support, patience, advice and guidance that he has given me through out
these years.
Boundless thanks to my parents, David and Rosa, and my sisters, Nancy and
Fabiola, for their unfailing love, sacrifice and support that helps me endure even the
most difficult times. Thanks also to everyone in the UBF, my second family, for all
their love and selfless help that they have always offered me.
I would also like to thank my committee members: Dr. Panagiotis Dimi-
trakopoulos, Dr. Michael Zachariah, Dr. Jeffrey Klauda and Dr. Oded Rabin for
their insightful and valuable comments and suggestions.
I also want to thank the A. James Clark School of Engineering Fellowship for
funding my first year of the doctorate program, the US National Science Foundation
for providing funding for this project through grant CBET1160132, and the CONA-
CyT (National Council of Science and Technology of Mexico) for the scholarship
provided to complete my studies.
ii
Table of Contents
List of Figures v
List of Abbreviations ix
1 Introduction 1
1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.1 Copper and copper oxide thin films . . . . . . . . . . . . . . . 1
1.1.2 Modeling chemical vapor deposition (CVD) . . . . . . . . . . 4
1.2 Model reduction and stiffness relaxation background . . . . . . . . . . 5
2 Methodology 9
2.1 Reaction factorization on a lumped model deposition system . . . . . 9
2.1.1 Integration of the DAE system in time . . . . . . . . . . . . . 13
2.1.2 Outer solution of the singular perturbation problem . . . . . . 14
2.1.3 DAE index . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.1.4 High net-forward rates and multiple deposition reactions . . . 16
3 Reaction factorization on a spatially distributed deposition system 20
3.1 Simplified deposition system . . . . . . . . . . . . . . . . . . . . . . . 20
3.1.1 Tubular deposition reactor model . . . . . . . . . . . . . . . . 21
3.1.2 Reaction factorization and the singular perturbation problem . 24
3.1.3 Representative results . . . . . . . . . . . . . . . . . . . . . . 29
4 Application of the factorization approach to copper CVD 39
4.1 Copper chemical vapor deposition (Cu-CVD) . . . . . . . . . . . . . . 39
4.1.1 Surface reactions and growth model . . . . . . . . . . . . . . . 42
4.1.2 Gas phase reactions . . . . . . . . . . . . . . . . . . . . . . . . 44
4.1.3 Precursor depletion spatial profiles in the hot-wall tubular
CVD reactor . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.1.4 Application of the reaction factorization . . . . . . . . . . . . 48
4.1.5 Integration of the DAE system . . . . . . . . . . . . . . . . . 51
4.1.6 Simulation results and discussion . . . . . . . . . . . . . . . . 53
5 Numerical tools 59
5.1 Chemical species class . . . . . . . . . . . . . . . . . . . . . . . . . . 59
5.2 Adsorption isotherms from absolute rate theory . . . . . . . . . . . . 62
5.2.1 Copper adsorption isotherm . . . . . . . . . . . . . . . . . . . 63
6 Experimental deposition and image processing analysis of copper films 67
6.1 Experimental reactor system . . . . . . . . . . . . . . . . . . . . . . . 67
6.2 Results and analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
6.2.1 Experimental island size distribution . . . . . . . . . . . . . . 68
6.2.2 Average side length analysis . . . . . . . . . . . . . . . . . . . 70
iii
7 Conclusions and suggestions for future work 73
A Temperature and precursor influence in CuX O film growth and composition 80
A.1 Reactor system and substrate preparation . . . . . . . . . . . . . . . 80
A.1.1 Substrate preparation . . . . . . . . . . . . . . . . . . . . . . 81
A.1.2 Reactor operating procedure . . . . . . . . . . . . . . . . . . . 81
A.2 Representative films and design of experiments . . . . . . . . . . . . . 82
A.2.1 Representative results . . . . . . . . . . . . . . . . . . . . . . 83
A.2.2 Design of experiments . . . . . . . . . . . . . . . . . . . . . . 84
A.2.3 Image analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
A.3 Response surface modeling . . . . . . . . . . . . . . . . . . . . . . . . 88
A.3.1 Analysis of the full model . . . . . . . . . . . . . . . . . . . . 89
A.3.2 Reduced model . . . . . . . . . . . . . . . . . . . . . . . . . . 90
B Derivation of discretization weights 93
B.1 Second order accuracy Taylor’s expansion . . . . . . . . . . . . . . . . 93
B.1.1 Finite differences of first and second dervatives . . . . . . . . . 93
Bibliography 95
iv
List of Figures
1.1 Deposited films of Cu2O and CuO on two different types of substrates.
Figure A shows Cu2O (red color) and CuO (black color) deposited
on a 1 inch by 1 inch quartz substrate, where oxygen flow rate is
reduced from left to right. Figure B shows Cu2O deposited on a
copper substrate. Figure C shows a close up to the sharp transition
from Cu2O to CuO on the quartz substrates. The length of the white
line is 1 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.1 Schematic of the simplified deposition mechanism. Monomer precur-
sor M reversibly produces the trimer T in the gas phase, and can be
deposited through an irreversible reaction with a surface site X to
produce an adsorbed species A in the film. . . . . . . . . . . . . . . . 10
2.2 Two-dimensional dynamics with k0 = K0 = 1, σ = 1, and t ∈ [0, 2].
Specified initial conditions and those projected onto the equilibrium
manifold (green dashed curve Q) are marked with the filled blue square
and red circle, respectively. Each point denotes 0.02 time units; black
dotted lines denote [M]+3[T ] = constant. The blue curve corresponds
to  = 0.1 and the wto ODE model; the red trajectory represents the
outer solution of the singular perturbation problem. . . . . . . . . . . 16
3.1 Tubular reactor schematic. Precursor M evaporates from a precursor
boat (white rectangle) and flows toward the deposition zone 0 ≤
ζ ≤ 1. co represents the specified boundary conditions in which we
assume a constant concentration of M in the gas phase upstream from
the deposition zone. c0 is the projection of the specified boundary
conditions onto the manifold Q (defined in (3.5). c0
+
and c1 are the
concentrations resulting from boundary conditions at ζ0 and at ζ1,
respectively, in the two-point boundary value problem resulting from
the factorization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.2 Schematic of the collocation points in the ζ dimension. . . . . . . . . 31
3.3 Upper plot: Residuals during the spatial integration of the DAE sys-
tem (3.4) showing a superlinear convergence order of -1.626 (the norm
of residual axis is in logarithmic scale). Lower plot: Grid convergence
analysis of the solution to the system in (3.4) (the error axis is in log-
arithmic scale) with an order of convergence of -2.635. . . . . . . . . . 36
v
3.4 Phase space plot of the system. The equilibrium manifold Q is shown
as the solid blue curve and the spatial evolution of the system is
graphed as the light red section along the manifold; each point in
this section represents the composition in the reactor at a different
point in space ζ. The boundary conditions are indicated by the yellow
points. For reference on these conditions see Figure 3.1. The green
plane represents constant x0, i.e. precursor conservation and so any
co in this plane will be projected to c0 as a result of conservation of
species. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.5 Profiles of the concentration of monomer, [M], trimer, [T], and ad-
sorbed, [A], species along the reactor length ζ. . . . . . . . . . . . . . 38
3.6 Integrated profile of the surface property L, film thickness in this
example, along the reactor coordinate ζ at different τ times. . . . . . 38
4.1 Schematic of the hot wall tubular CVD reactor. The inner tube is
not used for this deposition process. Argon flows through the outer
tube. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
4.2 Image analysis process of copper films deposited on quartz substrates
using CuI as as copper precursor. The image on the left is a picture
of the copper (red) deposited on a 1 inch by 1 inch quartz substrate.
The image at the center is a optical microscope image of an area of
the film. The image on the right is the digitally process image of the
optical microscope image at the center. The processing was realized
by filtering the microscope image digitally through a threshold value
of 127 in the RGB color scale (the middle of the scale from 0(black)
to 255(white) and converting it to black and white. Note the accumu-
lation of copper clusters on the scratches of the film where the larger
clusters form, a potential indicator of a surface reaction deposition
mechanism. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
4.3 Sketch of the proposed growth mechanism. Gold color cubes rep-
resent copper atoms and red cubes iodine atoms. Copper can be
incorporated to a growing island through direct impingement of pre-
cursor molecules, Cu and CuI, from the gas phase, or through surface
diffusion of copper adatoms (yellow color). . . . . . . . . . . . . . . . 41
4.4 Sketch of the top view of the area around the edge of the growing Cu
crystal (in yellow) on a quartz surface (in light blue). . . . . . . . . . 43
4.5 The gas-phase reactions and their connection to the deposition pro-
cess. D(s) is a deposition site that can be either a copper surface site
or a quartz surface site. All species are gas phase species unless noted
as surface (s) species. . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.6 Plug flow sketch for the precursor depletion model. . . . . . . . . . . 46
vi
4.7 Concentration profiles of precursors along the reactor. The spatial
integration of the system is done on the deposition zone. This zone is
the shaded area in the reaction schematic of Figure 4.1. The boundary
condition at the beginning of the deposition zone is indicated by the
filled circle in the copper iodide profile. The light blue shaded area
is the physical position of the substrate on the reactor coordinate. . . 56
4.8 Average island size in the growing film along the reactor dimension.
The light blue shaded area is the physical position of the substrate
in the reactor coordinate. Note that the model predicts the average
size length of the islands to match the experimental average island
size in the substrate. . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.9 A Snapshoot at time = 5 minutes of the deposition rates of copper
and copper iodide on the surface of growing islands. The plot on the
left compares direct deposition rates from the gas phase, and the plot
on the right compares deposition rates by surface diffusion, The scale
of the rate axes is logarithmic. . . . . . . . . . . . . . . . . . . . . . . 58
5.1 Representation of the species classes. atomSurf is a simple class that
represents a surface site, it requires only the molecular weight of the
species, its density, and it calculates the site density [Xˆ ]. . . . . . . . 60
5.2 Schematic of an adsorption reaction. A is the precursor molecule in
the gas phase, X is a surface site available for deposition and AX is
the adsorbed molecule A on the site X . . . . . . . . . . . . . . . . . 63
5.3 Fractional coverage, θ, of copper on deposited copper and on a quartz
substrate as a function of the partial pressure of copper in the gas
phase. The red vertical lines indicate the copper partial pressure at
the reactor temperature, PvaCu = 0.133× 10
−4 Pa, and the copper
partial pressure at which our reactor operates, PCu = 4 × 10 −4 Pa. . 66
6.1 Reproduction of Figure 4.2. The image on the left is a quartz sub-
strate with copper deposited(red color), the image on the center is a
microscope image of an zone in the substrate, and the image on the
right is the black and white version of the microscope image. Details
on the image analysis process are described in this section and also
in the caption of 4.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
6.2 Overall island size distribution. The size of the islands was estimated
taking the pixel area as the base area of a cubic island. . . . . . . . . 71
A.1 CuI CVD reactor process tube and furnace (left); schematic of CVD
reactor showing relative positions of CuI precursor boat, O2 mixture
feed tube and substrate location. L1 = 18cm and L2 = 32cm. . . . . . 81
vii
A.2 Cu2O/CuO deposition on quartz (A); O2 mixture flows for the three
samples are M = 60 (left), 30 (center), and 20 (right). Cu2O on Cu
substrate is shown as (B) corresponding to M = 30 O2 mixture flow.
The grid shown in the background of A and B corresponds to 1 cm
spacing. Image (C) illustrates the sharpness of the transition between
the films; the horizontal white bar corresponds to 1 mm. This film
was deposited with M = 15. For all cases, deposition was performed
at T = 750 oC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
A.3 The image processing approach to determining film composition from
original film images (top). Note the bimodal distribution in the his-
togram (center), indicating the fraction of each form of the copper
oxide. Pie charts at the bottom indicate relative ratios of Cu2O (red)
and CuO (black). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
A.4 Full model predictions for film composition (left) and the reduced
model (right) compared to modeling (blue circles) data; dashed red
lines indicate ±σM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
viii
List of Abbreviations
, δ Perturbation parameters
ζ Spatial dimension
t Time
Eg ,i Band gap energy of species i
G Net-forward reaction rate associated with a fast mode
F Net-forward reaction rate associated with a slow mode
[M] Concentration of precursor monomer
[T ] Concentration of precursor trimer
[X ] Concentration of available surface sites
[A] Surface concentration of adsorbed material A
[B] Concentration of bulk phase
L Film property whose evolution is tracked during deposition simulations
Kj Equilibrium constant for reaction j
kj Rate constant for reaction j
c Composition state vector
S Stoichiometric reaction coefficients matrix
r Net-forward reaction rates vector
nc Number of chemical species
nr Number of reaction rates
U Coordinate transformation matrix
z State vector of new coordinates associated with fast dynamic modes
x State vector of new coordinates associated with slow dynamic modes
w State vector of new coordinates associated with conserved quantities
y State vector of new reaction coordinates
Ux Row vectors of U corresponding to coordinates x
Uw Row vectors of U corresponding to coordinates w
g Vector with the equilibrium relationships of reversible reactions
f Vector of slow reaction rates
σ Surface area to volume ratio
Q Equilibrium manifold
∆t Size of time step used in Euler integration
JT Flux of trimer species in the gas phase
JM Flux of monomer species in the gas phase
DT Diffusion coefficient of the trimer species in the gas phase
DM Diffusion coefficient of the monomer species in the gas phase
v Average flow velocity
res Vector containing the residuals of the discretized equations in Chapter 3
E Discretization array with collocation weights for second order derivative
F Discretization array with collocation weights for first order derivative
F0 Discretization array with the weights for the first derivative at inlet
Fin Discretization array with the weights for the first derivative at interior points
FN Discretization array with the weights for the first derivative at outlet
Ein Discretization array with the weights for the second derivative at interior points
N Number of collocation points
ix
error Vector of error between solutions at different grid sizes
yfine Results form simulations with a fine grid in Chapter 3
yNi Results form simulations with a grid with Ni number of collocation points
τ = δt Slow time scale
kB Boltzmann’s constant
Avo Avogadro’s constant
h Planck’s constant
θ Fractional surface coverage
ηCu Sticking coefficient of Copper on quartz substrates
ηCuI Sticking coefficient of Copper Iodide on quartz substrates
rCug Rate of deposition by direct impingement of Cu from the gas phase
rCuIg Rate of deposition by direct impingement of CuI from the gas phase
rCus Rate of deposition by incorporation from the edge sites
rCuIs Rate of deposition by incorporation from the edge sites
Ag Copper island surface area available for deposition
Aedge Surface area available for deposition on the edge of a Cu island
LCu Length of the side of a solid copper island of one atom
nCu Number of atoms in a growing copper island
nCuside Number of atoms in a side of a growing copper island
L Average islands size length (in Chapter 4)
mCu Weight of an atom of copper
mCuI Weight of an atom of copper iodide
vCu1 Atomic volume of copper
vSiO2 Molecular volume of SiO2
ˆ[X ] Surface site density of quartz
Ro Inner radius of tubular reactor
A Area perpendicular to flow direction (Chapter 4)
Vζ Volume of differential element of analysis
∆ζ Size of spatial integration step
T Temperature of deposition experiments
nucSites Surface concentration of nucleation sites
[Cu3I3] Gas phase concentration of copper iodide trimer
[I2] Gas phase concentration of iodine
[Cu] Gas phase concentration of copper
[Cu2] Gas phase concentration of copper dimmer
[CuI] Gas phase concentration of copper iodide
ξel Energy level in electronic energy modes
ωel Degeneracy of energy level in electronic energy modes
ν Vibrational frequency energy for vibrational energy modes
σ Symmetry number for rotational energy modes (only in Chapter 5)
I’s Moments of inertia (Chapter 5)
ρ Density of bulk phase materials in growing films
Do Bond energy between adsorbed species and surface
trans() Species class method to calculate the translational partition function
x
elect() Species class method to calculate the electronic partition function
vibra() Species class method to calculate the vibrational partition function
rota() Species class method to calculate the rotational partition function
partfunc() Species class method to calculate partition function
pftv() Method to calculate partition function without the electronic contribution
A(T,P) Species method to estimate the Helmholtz free energy
U(T) Species method to estimate the internal energy
H(T) Species method to estimate the enthalpy
S(T,P) Species method to estimate the entropy
G(T,P) Species method to estimate the Gibbs free energy
monoatomG Species class for monoatomic species in the gas phase
polyatomG Species class for polyatomic species in the gas phase
atomSurf Species class for surface sites
adatom Species class for adsorbed atoms
admolec Species class for adsorbed molecules
ZA Partition function of gas-phase species A
ZX Partition function of surface site X
ZAX Partition function of adsorbed species A on site X
ν0 Vibration frequency of the bond between the adsorbed species and the surface
Cu−Cu Adsorption energy of copper on a copper surface (only in Chapter 5)
Cu−SiO2 Adsorption energy of copper on a quartz surface (only in Chapter 5)
∆3 Volume of the adsorbed species
PA Partial pressure of gas species A (in Chapter 5)
pCu Partial pressure of copper (in Chapter 5)
pvapCu Vapor pressure of copper (in Chapter 5)
K spads Adsorption isotherm constant of species sp
Vp Volume of a cubic copper island
Vm3 Volume of a cubic pixel
cpix2m Scaling factor to convert pixels to length
vCu Atomic volume of copper (in Chapter 6)
DepLen Length of deposition zone in the CVD reactor
ρCu Density of copper
L1 Length of the oxygen feed tube in the CVD reactor
L2 Length of the zone downstream from the oxygen feed in the CVD reactor
CuxO Copper oxide: cupric (x = 1), cuprous (x = 2)
DAE Differential-Algebraic System of Equations
ODE Ordinary Differential Equation
CVD Chemical Vapor Deposition
MOCVD Metal-Organic Chemical Vapor Deposition
ALD Atomic Layer Deposition
MW Molecular Weight
MWa Molecular Weight of adsorbed species
MWs Molecular Weight of surface species
CIGS Copper Indium Gallium Sulfide
xi
RGB Red, green, and blue color model
scipy Scientific python library
FCC Face Centered Cubic crystal structure
NIST National Institute of Standards and Technology
CCCBD Computational Chemistry Comparison and Benchmark Database
ITO Indium Tin Oxide
XRD X-Ray Diffraction
RSM Response Surface Model
kMC kinetic Monte Carlo
MATLAB Matrix Laboratory - Programming language developed by MathWorks
PYTHON Programming language created by Guido van Rossum at CWI
CWI Centrum Wiskunde & Informatica
Appendix A
M Vector with rotameter settings
m Scaled rotameter settings
t Scaled temperatures
dof Degrees of freedom
SSE Sum of squared residuals
SST Total sum of squares
x Vector with the scaled factors
nM Number of sets of experimental conditions
Y Percentage of Cu2O on the deposited film observed experimentally
yRSM Percentage of Cu2O estimated by the RSM
b0 Coefficient in the RSM determined from experimental data
b Vector of coefficients in the RSM determined from experimental data
B Matrix of coefficients in the RSM determined from experimental data
c Vector of coefficients of the RSM determined from experimental data
X Rectangular array with factor values for evaluation of the RSM
Y¯ Mean measured value of the modeling data set
R2 Coefficient of determination
σ2M Variance of the modeling error
p Number of identified parameters in the model
σ2b,k Parameter variance
s2b,k Approximation to the true parameter variance
tk Student’s t-test t-value k
α Confidence level for Student’s t-test
tinv MATLAB function to calculate the t-value of a given probability
Appendix B
fi Value of approximation of function f at collocation point i
ζi Value of the spatial coordinate at collocation point i
∆ζi+1 Spatial difference between collocation points i + 1 and i
xii
∆ζi−1 Spatial difference between collocation points i − 1 and i
ξi Denominator of collocation weights at collocation point i
N Total number of collocation points
xiii
Chapter 1
Introduction
This chapter presents the motivation to study the chemical vapor deposition
(CVD) of copper and copper oxide thin films due to their current and potential
applications as materials for solar energy harvesting, as well as the challenges en-
countered in modeling this process. A literature survey on model reduction and
stiffness relaxation is also provided in the second part of this chapter.
1.1 Motivation
1.1.1 Copper and copper oxide thin films
Copper is widely used as an interconnect material in electronic circuits [1, 2]
and is a key component used in absorbent layers in CIGS solar cells and photovoltaic
devices [3, 4]. It has been used in metal-halide discharge lamps in the form of copper
iodide [5, 6, 7], as well as a precursor, in CuI, in CVD of copper oxide thin films [8, 9].
The most common copper oxides, CuO and Cu2O, are natural p-type semi-
conductors that have shown potential for photovoltaic and photoelectrochemical
applications [10]. The abundance of copper and the practically null toxicity of the
oxides makes CuxO a very attractive material for solar energy applications. Par-
ticular interest has been devoted to solar water splitting applications because the
oxide bandgaps (Eg ,Cu2O ≈ 2.1 eV and Eg ,CuO ≈ 1.2 eV) are suitable for the photo-
1
splitting of water to produce hydrogen and oxygen using visible light [11]. Practical
application of CuxO to water-splitting must, however, overcome barriers such as
semiconductor film instability in electrolyte solutions either by oxidation of cuprous
oxide to cupric oxide [12] or by reduction of CuO to Cu2O [10].
In the majority of the applications mentioned above, solar cells and photoelec-
trochemical cells, the films and electrodes are created by chemical vapor deposition
(CVD) using a copper containing precursor such as [13] or by co-evaporation using
elemental copper in physical vapor deposition as in [4]. Previous studies of copper
oxide CVD using copper iodide as the Cu precursor [8, 9] have revealed unexpected
film deposition spatial patterns: both Cu2O and CuO deposit are found on the
same substrate separated by a sharp transition between the two oxides. The abrupt
change in film composition is clearly seen due to the red color of Cu2O and black
color of CuO as shown in Figure 1.1.
In earlier work [14], the temperature and oxygen partial pressure influence on
these spatial patterns were statistically analyzed by creating a design of experiments
with three levels in each variable and then identifying surface response models from
the experimental data; more details of this study can be found in the Appendix A.
Subsequent analysis of thermodynamic degrees of freedom for a system consisting of
both oxides, Cu2O and CuO, revealed that these patterns are due to the coexistence
of both copper oxides at a unique critical oxygen partial pressure for a given tem-
perature. To complete this study, a one-dimensional oxygen transport model was
2
A BC
Figure 1.1: Deposited films of Cu2O and CuO on two different types of substrates. Figure A shows
Cu2O (red color) and CuO (black color) deposited on a 1 inch by 1 inch quartz substrate, where
oxygen flow rate is reduced from left to right. Figure B shows Cu2O deposited on a copper substrate.
Figure C shows a close up to the sharp transition from Cu2O to CuO on the quartz substrates. The
length of the white line is 1 mm.
created to estimate the location of the sharp Cu2O-CuO transition on the films. The
model was able to predict qualitatively the influence of the temperature and oxygen
composition on the spatial patterns observed during the design of experiments.
The next step in modeling this deposition system was to include a deposition
mechanism for copper and oxygen on the substrates in combination with the species
transport equations. This not only increased the overall complexity of the deposition
model, but introduced the additional requirement of computing kinetic constants
unavailable in the literature. Furthermore, while operation of this CVD reactor may
appear to be at steady state, the film is growing under non-equilibrium conditions.
Therefore the need to simulate this system with at least two different time scales
combined with the lack of a some kinetic information motivated the creation of the
methodology presented in this thesis.
3
1.1.2 Modeling chemical vapor deposition (CVD)
Chemical vapor deposition (CVD) processes are widely used in the semicon-
ductor manufacturing industry. In this process, a gaseous mixture of precursor
species containing the desired materials for the film flows over a substrate where the
precursors deposit and react with the surface to form a thin film.
A detailed kinetic mechanism is vital in modeling the complex chemical re-
action system of CVD processes for design, optimization and control. However,
a detailed kinetic model involves a large number of kinetic parameters to be de-
termined, either from literature or experimentally. Furthermore, the scales of the
kinetic parameters can inherently differ from each other over several orders of mag-
nitude making the model stiff and computationally demanding.
The computational cost increases with the number of precursor species, cre-
ating a high-dimensional system, especially when gas phase and surface transport
phenomena are coupled to the gas and surface reacting rate expressions. It is there-
fore advantageous to develop a framework to reduce the order of these chemical
reaction systems models, while enabling the tracking of the evolution of the compo-
sition of species along the reactor in a time scale that allows sufficient accuracy for
prediction and control.
In this thesis, a systematic model reduction methodology that is independent
of the form of the reaction rate models is proposed. The methodology is based
on a reaction Gauss-Jordan factorization procedure that allows a recasting of the
modeling equations into a singular perturbation problem. The methodology is tested
4
first with simple deposition models and then is applied to the copper chemical vapor
deposition processes.
1.2 Model reduction and stiffness relaxation background
Some of the first researchers to describe a methodology for reducing models
of reaction systems were Mass and Pope in [15] and Lam and Goussis in [16]. The
objective of Mass and Pope was to simplify complex chemical kinetics by reducing
the dimensionality of the state space; Lam and Goussis aimed to decouple slow and
fast time scales and overcome the stiffness of the differential equation system. While
their goals in model reduction were different, they both used a similar approach
based on the concept of dynamical systems.
On the one hand, Mass and Pope noted that when solving a complete model
with detailed kinetics the trajectories in the composition state space tend to ap-
proach one another before equilibrium is reached creating an attracting manifold.
From this observation they developed a perturbation analysis of the eigenvalues of
the Jacobian of the system of differential equations to identify the trajectories for
which certain eigenvectors vanished. These vanishing trajectories were associated
with fast time scales that were then eliminated to reduce the model to only the
variables that define a manifold. The reaction system then evolves on this manifold
that roughly represents states that are in pseudo-equilibrium with respect to the
fastest relaxing time scales.
On the other hand, Lam and Goussis used a computational singular pertur-
5
bation approach to find a refined a basis vector to separate the time scales. The
basis vector is defined so that the vector product of the stoichiometric vector and
the reaction rate vector is diagonal or block diagonal and thus can separate the time
scales of the system. The method developed in [16] was later extended to applica-
tion in reaction-diffusion systems in [17] to try to reduce stiffness due to reaction
and diffusion processes. In the methodology proposed in this thesis, instead of an
eigenvector analysis of the Jacobian of the system or refining a basis vector, the
separation of times scales is carried out with a Gauss-Jordan factorization.
As seen in the works [15] and [16], the singular perturbation theory provides
a natural framework for model reduction. In general if the system can be mod-
eled in the standard singularly perturbed form, where the time scales are explicitly
separated by a small parameter , then the system can be reduced to a system of
differential-algebraic equations (DAE). However there is a wide range of systems
that can exhibit time-scale multiplicity in which the slow and fast dynamics cannot
be associated with specific process variables. In these cases it is necessary to use
a coordinate change that would allow explicit separation of the time scales. Re-
searchers Kumar, Christofides and Daoutidis proposed such a coordinate change to
model chemical processes with two time scales in [18].
In their proposed methodology, Kumar et. al. define the algebraic variables
as the limit of the physical constraints in the slow time scale at the limit → 0, and
base the coordinate change on these algebraic variables and the Lie derivative of the
physical constraints. Furthermore, they provide necessary and sufficient conditions
for the existence of a -dependent and of an -independent coordinate change based
6
on the involutivity1 of the matrix that weights the algebraic variables.
Vora and Dautidis [19] provided further insight into the construction of the
coordinate changes for nonlinear systems. To do this they defined a generalized
stoichiometric matrix that contained the stoichiometric coefficients of the chemical
reactions in the system and separate this matrix into two matrices containing the
columns for fast and slow reactions, respectively. Then they enforced linear inde-
pendence on the matrix with the stoichiometric coefficients of the fast reactions as
well as on the vector containing the reaction rates for the fast reactions. Conse-
quently they noted that the involutivity condition of the distribution spanned by
the columns of the stoichiometric matrices for the fast reactions is readily satisfied
for isothermal systems since the matrix does not contain state variables. Therefore
the new coordinates were a linear combination of the original coordinates where the
transformation coefficients belong to the null space of the transpose of the matrix
with the fast reactions stoichiometric coefficients. They also noted that the algebraic
equations in the DAE representation of the problem obtained after the coordinate
transformation corresponded to reaction equilibrium or complete conversion con-
straints.
Contou-Carrere and Daoutidis developed a model reduction method for reaction-
convection processes with the goal of eliminating the stiffness from the systems
in [20]. They used a combination of the method of characteristics and singular
perturbation theory because modal decomposition techniques cannot be applied to
1Involutivity is the property of a matrix to be its own inverse. e.g A2 = I where I is the identity
matrix
7
reaction-convection systems since all the eigenmodes of the spatial differential oper-
ator in these systems have the same amount of energy. After applying the method
of characteristics, they use singular perturbation theory to separate spatial from
temporal variables.
In the methodology presented in this thesis the coordinate change is provided
through the reaction factorization. This factorization, although simple in nature,
provides a coordinate transformation matrix that identifies conserved quantities
directly. Likewise, the approach developed in this thesis provides a rigorous means
of separating modeling time scales, resulting in a “nested” sequence of singular
perturbation problems for transient deposition models that are spatially dependent.
These and additional benefits and features of the reaction factorization method will
be described in the following chapters of this thesis.
8
Chapter 2
Methodology
This chapters describes the reaction factorization methodology for a lumped
deposition model. The ability of the factorization to guide the separation of times
scales is explored through a simplified deposition reaction system and a procedure
to solve the singular perturbation problem resulting from the reaction factorization
is outlined.
2.1 Reaction factorization on a lumped model deposition system
To illustrate the methodology proposed in this thesis, consider a simplified
deposition mechanism consisting of a gas-phase monomer M and a gas phase trimer
T species and the reversible reaction between the two1. A molecule of precursor
monomer M, deposits a single atom of material A through an irreversible reaction
with a surface site X producing the growth of the film, as sketched in Figure 2.1
where the net-forward reversible and deposition reaction rates are represented by
G0 and F0, respectively. In this analysis it is assumed that the net-forward rate
associated with the reversible reaction is much greater than the rate of the deposition
reaction. Secondary products of the surface reaction are not considered in this
1While the choice of a trimer species may seem at first unusual relative to a dimer, we will show
later that Cu3I3 is the only oligomer formed in the gas phase from CuI
9
M TMM+ +
A X X X XA A X A XX X A X A 3M G0M + X TAF0
G0 =
1

(
K0 [M]
3 − [T]
)
mol m−3 s−1
F0 = k0 [M] [X] mol m−2 s−1
Figure 2.1: Schematic of the simplified deposition mechanism. Monomer precursor M reversibly
produces the trimer T in the gas phase, and can be deposited through an irreversible reaction with
a surface site X to produce an adsorbed species A in the film.
simplified analysis. Material balances for each species yield the following ordinary
differential equations in time:
dcnc×1
dt
= Snc×nr rnr×1 =⇒
d
dt












[M]
[T ]
σ [X ]
σ [A]












=












−3 −1
1 0
0 −1
0 1
















G0
σF0



 (2.1)
where σ is the deposition surface area to reactant gas volume ratio with units m−1.
The matrices c, S, and r contain the species concentrations, stoichiometric coeffi-
cients, and reaction rates respectively. nc and nr refer to the number of chemical
species and reactions. The species concentrations are given in mol/m3 for gas phase,
[M] and [T], and mol/m2 for [X] and [A].
This original formulation (2.1) of the problem presents a higher-dimensional
set of ODEs, four equations, compared to the single rate-limiting deposition step.
Moreover, the dynamic behavior of the monomer concentration, [M] (t) is governed
by two reaction rates of different orders of magnitude.
It might initially appear that the QR or singular value decomposition methods
10
would be the best candidates to diagonalize the system of differential equations 2.1
and to subsequently reduce it to its minimal dynamic dimension. Nevertheless, these
methods are not useful for guiding the separation of time scales because the matrix S
is defined by reaction stoichiometric coefficients instead of reaction rate coefficients.
On the other hand, a Gauss-Jordan factorization of S can provide guidance in the
time scale separation by pairing each reaction with each new reaction coordinate to
the greatest extent possible, similar to the approach of Vora and Daoutidies in [19].
Hence, dynamic dimension reduction and time scale separation in Equation 2.1 can
be obtained if a transformation matrix U exists such that in the coordinate change
y = Uc each new coordinate can be associated to each reaction rate, yi ∈ Rnc , i.e.,
Unc×ncSnc×nr ≈




Inr×nr
0(nc−nr )×nr



 .
Clearly, complete diagonalization of the stoichiometric matrix should not be ex-
pected in every case. However, this factorization can still guide the time scale
separation even if off-diagonal elements remain. For the system under study in this
section, the ODEs in Equation 2.1 can be factored to find
dy
dt
=
d
dt












[T ]
− [M]− 3 [T ]
σ [X ]− [M]− 3 [T ]
σ [A] + [M] + 3 [T ]












=
d
dt












z0
x0
w0
w1












=












1 0
0 1
0 0
0 0
















G0
σF0



 (2.2)
yielding the following temporal modes for the specific case k0 = 1 m3mol−1s−1,
11
K0 = 1 m3mol−1, and σ = 1 m−1:
fast (reversible)
dz0
dt
=
d [T ]
dt
=
1

(
[M]3 − [T ]
)
(2.3)
slow (deposition)
dx0
dt
=
d(−[M]− 3[T ])
dt
= [M][X ] (2.4)
conserved quantities w0 = [X ]− [M]− 3[T ] (2.5)
w1 = [A] + [M] + 3[T ] (2.6)
or more generally, when the factorization 2.2 can be completed,
slow
dx
dt
=
dUxc
dt
= F(c)
dx
dt
= F(x, z,wo) (2.7)
fast
dz
dt
=
dUzc
dt
=
1

g(c) →0=⇒ 0 = g(x, z,wo) (2.8)
conservation w = Uwc (2.9)
where G = g/. We observe a semi-explicit differential-algebraic (DAE) system
results, both in the form of example (2.3-2.6) and in the more general case (2.7-2.9);
algebraic equations result from both the elimination of reduntant dynamic modes
producing the conserved quantities, and as will be discussed later, as → 0.
Following the nomenclature by Biegler in [21], x are the differential variables,
z the algebraic variables, and yT = [xT zT wT ] [21]. Ux ,Uz ,Uw are submatrices
from U containing the rows corresponding to the slow, fast and conservation new
coordinates, x, y, z, respectively. It should be observed that new coordinates w
reveal quantities conserved during the time evolution the slow modes. For example,
12
w1 is simply the total number of deposition element atoms in our control volume
and w0 + w1 = [X ] + [A] corresponds to the conservation of surface states. It is
important to recall that the decomposition approach presented here resembles other
model reduction procedures such as the ones in [18, 19, 22, 23] with the differences
noted in the introduction section of the thesis.
2.1.1 Integration of the DAE system in time
The system of differential-algebraic equations resulting from the reaction fac-
torization in equations 2.4-2.6 can be integrated in a relatively direct manner for
a finite . The ODEs involving x and z are integrated in time by computing the
reaction rates in the original coordinates by first recovering c with the inverse of the
transformation matrix from c = U−1y. It is evident that U must be invertible since
the forward elimination procedure used to obtain it is reversible. It should also be
noted that the new coordinates representing the conserved quantities w = wo are
clearly constant for all time.
Figure 2.2 shows representative results for the DAE system in equations 2.4-
2.6 calculated using a fixed step-size implicit Euler integrator. The initial specified
conditions used are co = [0.4, 0.6, 1, 0]T resulting in xo = −2.2, zo = 0.6, and wo =
[−1.2, 2.2]T after the transformation. The figure illustrates the strong attraction
of the initial conditions co to the neighborhood of the equilibrium manifold Q and
the subsequently slower dynamics which follow. The nature and computation of
manifold Q are described next.
13
2.1.2 Outer solution of the singular perturbation problem
The reaction system takes the form of a singularly perturbed ordinary differen-
tial equation when Equation 2.3 is multiplied by . Furthermore, in the limit → 0
the fast dynamics in the species composition space vanish, relaxing the system to
the manifold defined by the equilibrium condition
Q =
{
[M], [T ] : K0[M]3 − [T ] = 0
}
. (2.10)
The numerical outer solution of the singular perturbation problem [24] with initial
conditions co is calculated in the following manner.
First, the initial specified conditions co = c(t = 0) are converted to the new
chemical species composition space using the transformation matrix U, i.e yo = Uco .
Then these initial conditions co are projected onto the equilibrium manifold Q
defined by (2.10), to find c0 by solving the set of nonlinear/linear algebraic equations:








Uxc0
g(c0)
Uwc0








=








xo
0
wo








.
It is important to notice that while the manifold Q is defined in the original coordi-
nates, i.e. [M] and [T ], the above set of nonlinear/linear equations gives the initial
conditions in the new coordinate space y. The current state of the new coordinates
for the species state ci at time t i (where i = 0 for the projected initial state) is then
computed during the integration process using the transformation matrix Ux .
xi = Uxci .
14
Finally we use a Newton-Raphson or another suitable nonlinear algebraic equation
solver to calculate the implicit Euler update over time step ∆t by solving driving
the residual below to zero.








(xi+1 − xi )/∆t − f
(
U−1yi+1
)
g
(
U−1yi+1
)
wi+1 −wo








.
Figure 2.2 shows the results of applying this method to the system in equations
2.3-2.6 for  = 0. The figure shows the specified initial conditions co , the projection
of these conditions into the equilibrium manifold Q to give the initial conditions c0,
and the evolution in time of the system along the manifold approaching equilibrium
at the point c∞.
2.1.3 DAE index
It should be mentioned that when solving semi-explicit DAE systems one must
be concerned about the index of the system and the consistency of initial conditions.
The index of a DAE is defined as the number of times that the system must be
differentiated with respect to time in order to convert it into an explicit set of first
order ODEs [25]. For DAE systems with index higher than one, the differentiation
to reduce the index to one generates more algebraic variables. Consistent initial
values are those who satisfy not only the original algebraic variables but also the
ones resulting from the index reduction [26].
The semi explicit DAE system in equations 2.3-2.6 is index-1, thus consistency
of the initial conditions is not an issue in this case. Nonetheless, more complex
15
deposition surface reaction models may have higher DAE index values and thus
need more sophisticated treatment of the initial conditions [21, 27].
0.4 0.5 0.6 0.7 0.8
[M]
0.2
0.3
0.4
0.5
0.6
[T]
co
c0
Q
c
Figure 2.2: Two-dimensional dynamics with k0 = K0 = 1, σ = 1, and t ∈ [0, 2]. Specified initial
conditions and those projected onto the equilibrium manifold (green dashed curve Q) are marked
with the filled blue square and red circle, respectively. Each point denotes 0.02 time units; black
dotted lines denote [M] + 3[T ] = constant. The blue curve corresponds to  = 0.1 and the wto
ODE model; the red trajectory represents the outer solution of the singular perturbation problem.
2.1.4 High net-forward rates and multiple deposition reactions
Another numerical challenge encountered when solving models of surface re-
actions arises when the equilibrium reactions favor the product species and the de-
position reaction rate is significant. Under these conditions the numerical implicit
16
integrator may converge to negative values of reactant concentrations and provide
results with non-physical meaning. To illustrate this consider the following reaction
system
M + 2N
G0⇀↽ P
D + E
F0→B
with the net-forward reaction rates
d [M]
dt = −G0
d [N]
dt = F0 − G0
d [P]
dt = G0
where the reactants D and E are fed in such a way that the rate F0 is constant.
The matrix representation of the system is given by
d
dt








[P]
[N]
[M]








=








1 0
−2 1
−1 0












G0
F0




with G0 =
1

g0; after applying the factorization the system becomes
[M] + [P] = w0 the conserved mode
[N] + [P] = x0 + ∆tF0 the slow mode evolving over time interval ∆t
then the equilibrium manifold is defined by g0=0 and the integration of the slow
mode to give the equations:
[P]− K0[M][N] = 0
x0 + ∆tF0 − [N]− K0 (w0 − x0 −∆tF0 + [N]) [N] = 0
17
Setting the equilibrium constant K0 = 1 and solving for [N] results in
[N] =
x0 + ∆tF0 − 1− w0 ±
√
(1 + w0 − x0 −∆tF0)2 + 4(x0 + ∆tF0)
2
thus, even though only the positive solution for [N] has a physical meaning, if
K0  1 the numerical integrator can overshoot the limit [M] → 0 and converge to
the negative solution as mentioned above.
Finally, the presence of multiple, relatively slow deposition reactions that all
deposit the same chemical species may prevent the reaction factorization from fully
diagonalizing the stoichiometric matrix. However, as mentioned before, the result
may still be valid and useful for guiding the dimension reduction and separation of
time scales. For example consider adding the following deposition reaction to the
system of equations 2.1.
T + X
F1→ A + 2M F1 = k1[T ][X ] (2.11)
using the same values for the reversible and deposition reaction constants as in the
analysis above, k0 = K0 = 1, and k1 = 0 the resulting DAE is
slow
dx0
dt
=
d([M] + 3[T ])
dt
= −[M][X ]− [T ][X ]
fast
dz0
dt
=
d [T ]
dt
=
1

(
[M]3 − [T ]
)
− [T ][X ] →0=⇒ [M]3 − [T ] = 0
conserved w0 = [X ] + [A]
w1 = [A] + [M] + 3[T ].
In this case it can be seen that the conserved quantities are the same as those ob-
tained before the addition of the new deposition reaction (2.11): the conservation of
18
surface states and the conservation of the deposited material found in species M, T,
and A. The equilibrium relationship resulting is also the same, and a single dynamic
mode now describes the consumption of the total precursor by the combination of
both surface reactions. The kinetic model considered for copper CVD in Chapter
4 presents this case where the factorization does not fully diagonalize the system.
However, the model is structurally correct and the singular perturbation problem
will be shown to be valid.
19
Chapter 3
Reaction factorization on a spatially distributed deposition system
This chapter studies the reaction factorization of a spatially distributed de-
position system, extending the analysis in the previous chapter by considering gas
phase species diffusion in the hot wall tubular reactor. This model provides a more
realistic CVD process relative to the lumped parameter model in Chapter 2, and is
more closely related to our copper deposition system.
3.1 Simplified deposition system
To demonstrate the next step in extending the reaction factorization proce-
dure, consider a model CVD system in which (metal) precursor monomer species
M in the gas phase adsorbs onto the growth surface to form the adsorbed species
A with net adsorption rate G1; monomer M also can undergo a trimerization reac-
tion to form T at a rate G0. The adsorbed species is incorporated into the growing
film B through an irreversible reaction at rate F0; this same rate governs the time-
evolution of a film property L, such as film thickness, by the rate δF0, where δ relates
20
the concentration of deposited species B to the film property.
3M ⇔ T G0 =
1
0
(
K0 [M]
3 − [T]
)
mol s−1m−3
M ⇔ A G1 =
1
1
(K1 [M]− [A]) mol s−1m−2
A → B F0 = k0 [A] mol s−1m−2
(3.1)
The reactions and their associated rate expressions are summarized in Equa-
tion 3.1, where K0 is the equilibrium constant associated with the trimerization
reaction with units mol−2m6, K1 is adsorption equilibrium coefficient with units of
m, and k0 is the forward rate constant of the irreversible reaction with units of s−1.
The small parameters i  1 represent the time-scale ratios (with units of s) of
the reversible reactions relative to the irreversible reaction, which is assumed to be
the rate-limiting step of this example. These small parameters mark the separation
of fast and slow modes in the singular perturbation formulation of the problem; a
physical interpretation of this system is that the deposition process is kinetically
limited.
3.1.1 Tubular deposition reactor model
The deposition reactions will take place in a tubular reactor system, as the
one depicted in Figure 3.1, and thus a heterogeneous plug flow reactor (PFR) model
with axial diffusion and convection of the gas-phase precursors is derived; the axial
coordinate is denoted ζ. For this reactor system, the start of the deposition zone is
ζ = 0 and, without loss of generality, we set the reactor outlet as ζ = 1 m. As in
21
the previous chapter, surface and gas-phase molar concentrations are indicted with
the bracket notation [·].
=0 =1
co c0c0+ c1
Figure 3.1: Tubular reactor schematic. Precursor M evaporates from a precursor boat (white
rectangle) and flows toward the deposition zone 0 ≤ ζ ≤ 1. co represents the specified boundary
conditions in which we assume a constant concentration of M in the gas phase upstream from
the deposition zone. c0 is the projection of the specified boundary conditions onto the manifold Q
(defined in (3.5). c0
+
and c1 are the concentrations resulting from boundary conditions at ζ0 and
at ζ1, respectively, in the two-point boundary value problem resulting from the factorization.
The modeling equations are summarized in Table 3.1, where
Ji = −Di
∂ [i]
∂ζ
+ v [i]
is the flux (diffusion and convection) of species i . Only binary diffusion is considered
due to the low concentration of precursors. The parameter σ is the area-to-volume
ratio to convert surface concentration into volume based concentrations. Thus, for
a tubular reactor system of radius R where material is deposited on the inner wall
surface, σ = 2/R . The small parameter δ  1 captures the ratio of the deposition
time scale to the species transport time scale. Because of these small parameters,
the rates of change in the surface states can differ in orders of magnitude, rendering
the system numerically stiff.
22
Governing equation ζ = 0, t ≥ 0 ζ = 1, t ≥ 0 ζ ∈ (0, 1), t = 0
∂[T ]
∂t
=
∂JT
∂ζ
+ G0 JT (0) = u [T]0
∂ [T]
∂ζ
= 0 [T] = 0
∂[M]
∂t
=
∂JM
∂ζ
− 3G0 − σG1 JM(0) = u [M]0
∂ [M]
∂ζ
= 0 [M] = 0
∂[A]
∂t
= G1 − F0 −−−− −−−− [A] = 0
∂L
∂t
= F1 −−−− −−−− L = 0
∂[B]
∂t
= F0 −−−− −−−− [B] = 0
Table 3.1: Modeling equations for the heterogeneous tubular CVD system together with boundary
and initial conditions.
23
3.1.2 Reaction factorization and the singular perturbation problem
To eliminate redundant dynamic modes and to provide a rational path to
reconfiguring the modeling equations so that the system dynamics as 0, 1, δ → 0
can be investigated, the Gauss-Jordan factorization procedure proposed in Chapter
2 is applied to allow separating the effect of each reaction to guarantee the correct
handling of the different time scales of the model. Rewriting the modeling equations
in matrix form
∂
∂t
















[T]
[M]
[A]
L
[B]
















=
















∂JT/∂ζ
∂JM/∂ζ
0
0
0
















+
















1 0 0 0
−3 −σ 0 0
0 1 −1 0
0 0 0 1
0 0 1 0




























G0
G1
F0
F1












or
∂c
∂t
=
∂J
∂ζ
+ S












G0
G1
F0
F1












.
To find the transformation matrix U the factorization is applied to the stoi-
chiometric matrix, S, of the kinetic mechanism of our system yielding
24
SF = U · S =
















1 0 0 0 0
3 1 0 0 0
3 1 σ 0 0
0 0 0 1 0
3 1 σ 0 σ
















·
















1 0 0 0
−3 −σ 0 0
0 1 −1 0
0 0 0 1
0 0 1 0
















=
















1 0 0 0
0 −σ 0 0
0 0 σ 0
0 0 0 1
0 0 0 0
















. (3.2)
Multiplying the original state vector by the transformation matrix yields a set
of new coordinates given by
y = U · c =
















1 0 0 0 0
3 1 0 0 0
3 1 σ 0 0
0 0 0 1 0
3 1 σ 0 σ
















·
















[T]
[M]
[A]
L
[B]
















=
















[T]
3 [T] + [M]
3 [T] + [M] + σ [A]
L
3 [T] + [M] + σ [A] + σ [B]
















=
















z0
z1
x0
x1
w0
















(3.3)
After applying the factorization to the boundary and initial conditions, our
model transforms into the set of equations shown in Table 3.2. The table shows that
the model of the system in the new coordinates y, which are the linear combinations
of the original coordinates given in Equation 3.3. Note that the model equations
for each of the new coordinates involve only one net forward reaction rate, which
allows us to separate unambiguously the different time scale modes of the system.
The fast modes (for 0 < i  1) are tracked by the new coordinates z0 and z1,
which are the concentration of monomer and the total concentration of precursor
on the gas phase, respectively. Multiplying both equations by their respective small
parameter, 0 and 1, and taking the limits as the parameters go to zero, the outer
25
Governing equation ζ = 0, t ≥ 0 ζ = 1, t ≥ 0 t = 0, 0 < ζ < 1
∂z0
∂t
=
∂JT
∂ζ
+ G0 Jz0(0) = z0u
∂z0
∂ζ
= 0 z0 = 0
∂z1
∂t
=
∂ (3JT + JM)
∂ζ
− σG1 Jz1(0) = uz1
∂z1
∂ζ
= 0 z1 = 0
∂x0
∂t
=
∂ (3JT + JM)
∂ζ
− σF0 Jx0(0) = ux0
∂x0
∂ζ
= 0 x0 = 0
∂x1
∂t
= F1 −−−− −−−− L = 0
∂ (w0 − δx1)
∂t
= 0 −−−− −−−− L = constant + δ [B]
Table 3.2: Model structure after factorization.
26
solutions of the first two singularly perturbed PDEs are found by
g0 = 0 and g1 = 0
where Gi = gi/i .
The new coordinate w0 is a conservation mode, because the corresponding row
of the transformed stoichiometric matrix SF of Equation 3.2 vanishes, meaning it
has no generation or consumption terms. Furthermore, combining the variables x1
and w0 yields a conservation mode given by the fifth equation in Table 3.2, that
indicates that the final conserved quantity gives
L = constant + δ[B] = wo + δ[B]
which means the constant is w0 and represents the initial thickness or surface prop-
erty value.
The slow modes of the system, variables x0 and x1, correspond to the total
amount of precursor that has not been incorporated to the bulk of the film and to
the thickness of the film, respectively. The system evolves along these slow modes
during the slow time [28, 20]. Therefore defining the slow time τ = δt, taking
F1 = δF0 and recalling that x1 = L, we have that
δ
∂x0
∂τ
=
∂
∂ζ
(3JT + JM)− σF0
∂L
∂τ
= F0
therefore, when δ → 0
0 =
∂
∂ζ
(3JT + JM)− σF0
∂L
∂τ
= F0
27
In summary the evolution of the system in the slow time scale, along the manifold
defined by the equilibrium relations g0 and g1, when i → 0 and δ → 0, can be
tracked by the following equations
g0 = 0 (3.4a)
g1 = 0 (3.4b)
d
dζ
(3JT + JM) = σF0 (3.4c)
dL
dτ
= F0 (3.4d)
Our system is now described by a mixture of differential and algebraic equa-
tions, a DAE system, plus boundary and initial conditions to be described next.
To solve this system it is necessary to provide initial and boundary conditions that
are consistent not only with the algebraic equations, but also with the differential
equation 3.4c, in a manner referred to as consistent initialization [29].
The initial condition is defined as an empty reactor before the start of the
deposition, shown in Equation 3.6, which fulfills trivially the equilibrium manifold
given by
Q = {c : g0 = 0 and g1 = 0} (3.5)
in the gas phase
[M] (ζ, 0) = 0
[D] (ζ, 0) = 0
[A] (ζ, 0) = 0
L(ζ, 0) = 0
0 < ζ < 1. (3.6)
28
The specified boundary condition is co = [M0, 0, 0], in which we assume a con-
stant gas phase concentration of monomer on the boat zone as shown in Figure 3.1;
furthermore Equation 3.2 is used to obtain the specified conditions in the new co-
ordinate system. To translate these conditions to the slow time scale we project
them onto the manifold Q, this is the point c0 in Figure 3.1. Subsequently, c0 is
used to define the boundary conditions for Equation 3.4c, assuming that the gas is
moving at a uniform speed u. Since the system is constrained to the manifold Q,
equations 3.4a and 3.4b must be fulfilled together with the boundary conditions at
ζ = 0; this point is identified as c0
+
in Figure 3.1. The process for obtaining these
boundary conditions is summarized in Table 3.3 and a phase space representation
is given in Figure 3.4.
3.1.3 Representative results
As a representative simulation, consider an initially empty reactor with a M0 =
2 s−1m−3mol feed of pure monomer. Kinetic parameters for Equation 3.1 and trans-
port parameters used in this simulation are given in Table 3.4.
To solve our model described by Equation 3.4, one only needs to obtain the
steady state spatial profiles of precursors and adsorbed species in the reactor using
equations 3.4a through 3.4c, and then use this profile to integrate in time our tracked
property L over the slow time scale τ using Equation 3.4d. As mentioned in the
previous section, the first step in the solution of our model is to project the specified
boundary condition, co , onto the spatial distributions manifold to which the system
29
co c0 c0
+
c1
Variables
[M]o = M0 [M]
0 [M]0
+
[M]1
[T]o = 0 [T]
0 [T]0
+
[T]1
[A]o = 0 [A]
0 [A]0
+
[A]1
x0o = M0 x
0
0 x
0+
0 x
1
0
Equations
G0 = 0 G0 = 0 G0 = 0
G1 = 0 G1 = 0 G1 = 0
x0 = constant vx00 = Jx0
(
0+
) ∂
∂ζ
(x0) = 0
Table 3.3: Boundary conditions as they are transformed from the specified conditions co to the
projection onto the manifold Q, c0, and finally to the boundary conditions at the start (ζ0) and
end (ζ1) of the deposition zone, c0
+
and c1, respectively.
is constrained by the equilibrium relationships g0 and g1 and the transport equa-
tion 3.4c in Equation 3.1 when 0, 1, δ → 0. Then these projected conditions are
used to calculate the boundary condition at the beginning of the deposition zone;
these conditions are shown in Figure 3.4.
The second step in the solution of our system is the solution of the two-point
boundary value problem given by Equation 3.4c, the spatial differential equation
that drives transport in the reactor, while making sure the equilibrium relationships
in equations 3.4a and 3.4b hold for every integration point to ensure that the system
stays on the manifold of spatial distributions Q. For this, a set of collocation points,
30
Parameter Value Units
v 0.03 m/s
K0 1 m3mol−1
K1 1.5 m
k0 0.4 s−1
σ 1 m2m−3
DM 1.0×10−1 m2s−1
DD 1.0×10−2 m2s−1
Table 3.4: Kinetic and transport parameters for representative simulation of equations 3.4
ζi with i = 0, 1, ..., N , is defined in Figure 3.2 where ζ0 = 0 and ζN = 1. We then com-
0 i-1 i i+1 N... ...
Figure 3.2: Schematic of the collocation points in the ζ dimension.
pute discretization arrays E and F using second order accurate finite differences to
calculate the discretization weights of the second and first derivatives, respectively.
Derivation of the finite difference equations are given in Appendix B. Equation 3.4c
31
is then discretized at every point to evaluate the residuals in Equation 3.7
res =
















g0 (c)
g1 (c)
3 (DTF0− v) · T + (DMF0− v) ·M− v (3T 0 + M0)
3 (DTEin− vFin) · T + (DMEin− vFin) ·M− k0A
3 (DTFN) · T + (DMFN) ·M
















(3.7)
where the third and fifth equations are the residuals for the boundary conditions
at the start (ζ = 0) and end (ζ = 1) of the deposition zone respectively. F0 is
the first vector row in array F where the weights are calculated using forward finite
differences. This vector is used to calculate the residual of the boundary condition
at ζ = 0 for c0
+
in Table 3.3 (see Figure 3.4) which is the third element in the
residuals vector of Equation 3.7. F0 is written as:
F0 =
[
F(0, 0) F(0, 1) F(0, 2) · · · 0
]
where
F(0, 0) =
(ζ1 − ζ0)
2 − (ζ2 − ζ0)
2
ξF0
F(0, 1) =
(ζ2 − ζ0)
2
ξF0
F(0, 2) = −
(ζ1 − ζ0)
2
ξF0
with
ξF0 = (ζ1 − ζ0)
2 − (ζ2 − ζ0) (ζ1 − ζ0)
2 .
The matrix Fin is a matrix with the weights of the discretization of the first deriva-
32
tive at each interior collocation point. A vector row i of Fin is given by
F(i , i − 1) =
(ζi+1 − ζi )
2
ξF
F(i , i) =
(ζi−1 − ζi )
2 − (ζi+1 − ζi )
2
ξF
F(i , i + 1) =
− (ζi−1 − ζi )
2
ξF
where
ξF = (ζi−1 − ζi ) (ζi+1 − ζi )
2 − (ζi+1 − ζi ) (ζi−1 − ζi )
2
Fin is the equivalent of Ein for the discretization weights of the second derivative
for the interior collocation points. The weights of Ein for i = 1, 2, · · · , N are as
follows
Eini =
[
0 · · · E(i , i − 1) E(i , i) E(i , i + 1) · · · 0
]
i = 1, · · · , N
where central finite differences were used to calculate the weights as follow
E(i , i − 1) =
(ζi+1 − ζi )
ξE
E(i , i) =
(ζi−1 − ζi )− (ζi+1 − ζi )
ξE
E(i , i + 1) =
− (ζi−1 − ζi )
ξE
and
ξE =
(ζi+1 − ζi ) (ζi−1 − ζi )
2 − (ζi−1 − ζi ) (ζi+1 − ζi )
2
2
.
The vector FN is the last vector row of the F array. We use this vector to calculate
the boundary condition at ζ = 1 (column for c1 in Table 3.3). The vector has the
form
FN =
[
0 · · · F(N , N − 2) F(N , N − 1) F(N , N)
]
33
for this vector the weights are calculated with a backward finite difference to obtain
F(N , N − 2) =
− (ζN−2 − ζN−1)
2
ξFN
F(N , N − 1) =
(ζN−3 − ζN−1)
2
ξFN
F(N , N) =
(ζN−2 − ζN−1)
2 − (ζN−3 − ζN−1)
2
ξFN
with
ξFN = (ζN−2 − ζN−1) (ζN−2 − ζN−1)
2 − (ζN−3 − ζN1) (ζN−2 − ζN−1)
2 .
M = [M0, M1, M2, ..., MN ]
T are the discrete values of [M] at each collocation point
in space; similarly T, and A are the discretized values of [T] and [A] at each collo-
cation point. Finally c = [M,T,A]T . The residual (3.7) function was solved using
a Newton-Raphson method with N = 200 collocation points with a tolerance of
tol < 10−12 for the norm of the residual. The resulting profiles are shown in Fig-
ure 3.5 where it can be observed that the adsorption of the monomer on the surface
dominates the deposition and there is little trimer produced along the reactor, as is
expected by the thermodynamic and kinetic values employed in the simulation.
The plot in the upper part of Figure 3.3 shows the an analysis of the residuals
for a solution with N = 200. The linear behavior in the logarithmic scales (the
norm of residuals axis is in logarithmic scale) indicates that our system exhibits
superlinear convergence with an order of converge of -1.626. To analyze converge
with respect to discretization level, we compare solutions with different number of
collocation points against an analytical solution, or if an analytical solution is not
available, we compare against a fine grid solution. The number of collocation points
34
for each grid is calculated as Ni = 2ni +1 where ni = 2, 3, 4, · · · to generate grids with
matching collocation points and allow direct comparison of the function evaluations
at the same spatial location. Since an analytical solution is not available, we consider
a grid with Ni = 211 + 1 = 2049 collocation points to be our fine grid. The error
for a solution with Ni collocation points was computed as the norm of the difference
vector error = yfine − yNi . The result of the analysis can be seen in the lower part of
Figure 3.3, where the linear behavior of the residuals in both logarithmic scale axes
allows to compute an order of converge as the slope of the curve to be -2.635.
Finally to investigate the time evolution of the surface property L(ζ, τ), in
this case film thickness, L (ζ, τ) is integrated at every collocation point in the slow
time scale using Equation 3.4d to obtain the profile in Figure 3.6, where it has been
integrated over 100 τ units (m3s/#). Figure 3.6 shows the film thickness along
the reactor spatial coordinate, ζ at different τ times. Since in the slow time the
concentration profiles for precursors and adsorbed species are at quasi equilibrium,
the time integration of L consist of solving a simple set of ordinary differential
equations with a constant value (in time, but of course not in space) of [A] at every
collocation point. It can observed, that as expected, the thickness will grow more
quickly near the entrance of the reactor where the concentration of the adsorbed
species A is higher.
35
10-6 10-5 10-4 10-3 10-2 10-1 100 101
norm of update at iteration K
10-7
10-6
10-5
10-4
10-3
10-2
10-1
100
101
no
rm
 of
 up
dat
e a
t it
era
tio
n k
+1
ρ=0.99
101 102 103
number of collocation points
10-9
10-8
10-7
10-6
10-5
10-4
10-3
10-2
er
ro
r
ρ = -3.07
Figure 3.3: Upper plot: Residuals during the spatial integration of the DAE system (3.4) showing a
superlinear convergence order of -1.626 (the norm of residual axis is in logarithmic scale). Lower
plot: Grid convergence analysis of the solution to the system in (3.4) (the error axis is in logarithmic
scale) with an order of convergence of -2.635.
36
[M]0.0 0.5
1.0 1.5
2.0
[T] 0.0
0.1
0.2
0.3
0.4
[A]
0.0
0.2
0.4
0.6
0.8
1.0
1.2
co
c0
c0+ ζ=0
c1
ζ=1
Q
Figure 3.4: Phase space plot of the system. The equilibrium manifold Q is shown as the solid blue
curve and the spatial evolution of the system is graphed as the light red section along the manifold;
each point in this section represents the composition in the reactor at a different point in space ζ.
The boundary conditions are indicated by the yellow points. For reference on these conditions see
Figure 3.1. The green plane represents constant x0, i.e. precursor conservation and so any co in
this plane will be projected to c0 as a result of conservation of species.
37
0.0 0.2 0.4 0.6 0.8 1.0
ζ
0.00
0.05
0.10
0.15
0.20
0.25
0.30
0.35
0.40
0.45
[m
ol
/m
3
]
M
T
A
0.0 0.2 0.4 0.6 0.8 1.0
ζ
0
2
4
6
8
10
12
14
L
30
50
70
τ = 100 m3 s / #
Figure 3.5: Profiles of the concentration of monomer, [M], trimer, [T], and adsorbed, [A], species
along the reactor length ζ.
0.0 0.2 0.4 0.6 0.8 1.0
ζ
0.00
0.05
0.10
0.15
0.20
0.25
0.30
0.35
0.40
0.45
[m
o
l/
m
3
]
M
T
A
0.0 0.2 0.4 0.6 0.8 1.0
ζ
0
2
4
6
8
10
12
14
L
30
50
70
τ = 100 m3 s / #
t
L=δk0[A]
Figure 3.6: Integrated profile of the surface property L, film thickness in this example, along the
reactor coordinate ζ at different τ times.
38
Chapter 4
Application of the factorization approach to copper CVD
This chapter describes the application of the reaction factorization method-
ology presented in previous chapters to the modeling of a hot wall chemical vapor
deposition reactor. The factorization is shown to be able to guide the time scale
separation and to reveal conserved quantities in the new reaction coordinates.
4.1 Copper chemical vapor deposition (Cu-CVD)
This section presents the application of the reaction factorization methodology
to modeling of the chemical vapor deposition of copper for the reaction
CuI(g)→ Cu(s) +
1
2
I2(g)
in the reactor shown in Figure 4.1. Experiments on copper CVD and visual analysis
of deposited films are presented in Appendix B. From the experimental and digital
3.92 cm
50 cm
Ar
750oCCuI SiO2 substrate
18 cm
Figure 4.1: Schematic of the hot wall tubular CVD reactor. The inner tube is not used for this
deposition process. Argon flows through the outer tube.
image analysis to be presented in Chapter 6, a single-island growth was inferred and
39
is explained as follows: first, Figure 4.2 shows that the film grows mostly through the
formation, expansion and coarsening of islands, exhibiting a Volmer-Weber growth
mode. Presumably, this occurs because the copper adatoms are more attracted to
themselves than to the quartz surface, as suggested by the adsorption energies of
copper on copper (1.25 eV [30]) and copper on quartz (0.1 eV [31]). Hence, it is
assumed that the islands grow through two main precursor incorporation mecha-
nisms: 1) by direct precursor impingement from the gas phase onto the growing Cu
crystals; and 2) by initial adsorption onto the quartz, followed by surface diffusion
of the adsorbed species to the growing Cu crystals.
Figure 4.2: Image analysis process of copper films deposited on quartz substrates using CuI as as
copper precursor. The image on the left is a picture of the copper (red) deposited on a 1 inch by 1
inch quartz substrate. The image at the center is a optical microscope image of an area of the film.
The image on the right is the digitally process image of the optical microscope image at the center.
The processing was realized by filtering the microscope image digitally through a threshold value of
127 in the RGB color scale (the middle of the scale from 0(black) to 255(white) and converting it
to black and white. Note the accumulation of copper clusters on the scratches of the film where the
larger clusters form, a potential indicator of a surface reaction deposition mechanism.
40
3.
SiO
2
3.92
 
3.
2
c
3.
m 2 
3.m
2
c
3.m
50A9m
r
50A9m
r
Figure 4.3: Sketch of the proposed growth mechanism. Gold color cubes represent copper atoms and
red cubes iodine atoms. Copper can be incorporated to a growing island through direct impingement
of precursor molecules, Cu and CuI, from the gas phase, or through surface diffusion of copper
adatoms (yellow color).
41
4.1.1 Surface reactions and growth model
Figure 4.3 shows a schematic representation of a proposed copper crystal
growth mechanism in which precursor incorporation occurs by two main processes:
direct impingement from the gas phase, and surface diffusion of precursors deposited
on the substrate surface. In the first process, rCug and r
CuI
g , copper and copper io-
dide from the gas phase deposit directly on the the growing islands, resulting in
the addition of copper to the islands and the production of iodine. In the analysis
that follows, it is assumed that iodine atoms adsorbed on the Cu surface combine
instantly to form molecular iodide that immediately desorbs. The rate of this de-
position mechanism by direct impingement from the gas phase can be estimated by
the wall collision rates for gas phase Cu and CuI atoms using the kinetic theory of
gases and assuming an ideal gas
rCug = ηCu
√
kBT
2pimCu
[Cu] Ag , rCuIg = ηCuI
√
kBT
2pimCuI
[CuI ] Ag (4.1)
where kB is the Boltzmann constant, ηi is the sticking coefficient of precursor i , mCu
is the mass of a copper atom, mCuI is the mass of a copper iodide molecule, [Cu] and
[CuI ] are the gasp-phase concentrations of copper and copper iodide, respectively,
and Ag is the surface area available for deposition. To estimate this area, we consider
the islands to grow in a cubic manner. Therefore, for an island with nCu atoms of
copper, the length of the side of the island is L = n1/3Cu LCu where LCu = 0.228 nm is
the length of a solid copper atom represented as a cube. Then the area available for
deposition on the island is given by the following equation
Ag = 5× (L)
2 (4.2)
42
In the second mechanism the copper atoms deposit and weakly bond to the
quartz surface. These adatoms diffuse on the surface until they encounter the edge
of a growing island and are adsorbed by it. To estimate the rate of this adsorption
process through surface diffusion, absolute rate theory [32] is used to obtain equa-
tions 4.3 and 4.4, representing the rates of adsorption of Cu and CuI on the quartz
surface
rCus =
kBT
h
θCu/SiO2[Xˆ ]K
Cu
adsAedge (4.3)
rCuIs =
kBT
h
θCuI/SiO2[Xˆ ]K
CuI
ads Aedge (4.4)
where h is Plank’s constant, θCu/SiO2 is the equilibrium fraction coverage of copper
on quartz at the deposition temperature, [Xˆ ] is the surface site density of copper
on quartz, and Aedge is the area of edge sites available on the surface. The area
Edge sites
Figure 4.4: Sketch of the top view of the area around the edge of the growing
Cu crystal (in yellow) on a quartz surface (in light blue).
of edge sites, shown in Figure 4.4, available for deposition is estimated as follows:
first the number of copper sites in a side of the base of the island is calculated as
nCuside = (L/LCu + 1). Second the total area available for the sites on the perimeter
of the base is estimated as Aperi = nCusidev
2/3
Cu1 , where Cu1 is the volume of a solid
43
atom of copper. Finally to estimate the equivalent area of SiO2 surface sites, the
calculated copper edge area is divided by the area of a quartz surface site.
Aedge =
4
(
L
LCu
+ 1
)
Xˆ
v 2/3Cu1
v 2/3SiO2
To track the growth of the islands, we use the model given in Equation 4.5.
dL
dt
= δ
(
rCug + r
Cu
s + r
CuI
g + r
CuI
s
)
(4.5)
This completes the surface adsorption and growth model components.
4.1.2 Gas phase reactions
To calculate the concentration of precursors in the gas over the substrate
surface the following reactions are considered: (1) the reversible reaction between
the monomer CuI and the Cu3I3 trimer and (2) the reversible reaction between gas
phase atomic Cu and the Cu2.
The complete mechanism is shown in Figure 4.5, where D(s) represents an
surface adsorption site (either copper or quartz) and D(s)−Cu is an adsorbed copper
atom. The reactions and their rate expressions are presented in Table 4.1. In the
table, lower case k represents a forward reaction constant and upper case K an
equilibrium constant. The calculations of the adsorption and equilibrium constants
are detailed in Chapter 5.
44
Cu3I3 CuI I2 + 2Cu
CuD(s) + 0.5I2
Cu D(s)
Cu2
Figure 4.5: The gas-phase reactions and their connection to the de-
position process. D(s) is a deposition site that can be either a copper
surface site or a quartz surface site. All species are gas phase species
unless noted as surface (s) species.
0) 3CuI
k0K0⇀↽
k0
Cu3I3 r0 =
1
0
(
K0 [CuI]
3 − [Cu3I3]
)
1) 2CuI
k1K1⇀↽
k1
Cu2 + I2 r1 =
1
1
(
K1 [CuI]
2 − [Cu2] [I2]
)
2) Cu2
k2K2⇀↽
k2
2Cu r2 =
1
2
(
K2 [Cu2]− [Cu]
2
)
3) Cu
r3→ D(s) r3 = rCug ([Cu]) + r
Cu
s ([Cu])
4) CuI
r4→ D(s) + 12I2 r4 = r
CuI
g ([CuI ]) + r
CuI
s ([CuI ])
Table 4.1: The mechanism considered for Cu deposition from CuI. Reactions
1 through 3 are assumed to have relatively fast dynamics, where as reactions 4
and 5 the deposition reactions are considered to be the rate-limiting steps, and
i = 1/ki for i = 0, 1, 2.
45
4.1.3 Precursor depletion spatial profiles in the hot-wall tubular CVD
reactor
To compute the depletion behavior of the copper-containing precursors dur-
ing deposition, consider the reactor schematic in Figure 4.1 where the shaded area
represents where deposition is observed during experiments. The analysis in this
section assumes the gas flowing in the reactor to follow the plug flow model with
no diffusion. It is also assumed that there are no gas-phase composition gradients
in the radial direction. To derive the governing equations for this model, consider
the sketch in Figure 4.6 that models the interior of the reactor in Figure 4.1 where
E d
ge

Figure 4.6: Plug flow sketch for the precursor depletion model.
A = piR2o is the area normal to the flow and Ro = 1.95× 10
7nm is the inner radius of
the reactor tube, ci is the precursor species i concentration, v = 3×107nm/s (result-
ing from 570 sccm of Ar flowing at 750oC ) is the fluid velocity, and Vζ = piR2o ∆ζ
is the volume of the differential element of Figure 4.6. A mass balance over the
differential volume yields the following equation
piR2o ∆ζ
∂ci
∂t
= piR2o ci |ζν − piR
2
o ci |ζ+∆ζν + piR
2
o ∆ζri (c) i = 0, 1, .., 4
where ri is the overall rate of consumption of precursor i according to the rate ex-
pressions of Table 4.1 and ν is the average flow velocity. Dividing by the differential
46
volume and taking the limit when ∆ζ → 0 we obtain the deposition model
∂c
∂t
= −ν
∂c
∂ζ
+ r (c) (4.6)
where c and r are the vectors of concentration of species and reaction rate expres-
sions, respectively, as shown below
c =
















[Cu3I3]
[I2]
[Cu]
[CuI]
[Cu2]
















r =
















r0
r1 + 0.5r4
2r1 − 2r2 − σr3
−3r0 − 2r1 − σr4
r2
















.
Note that for the deposition reactions, r3 and r4 in Table 4.1, the rate expressions
must to be multiplied by the following surface area to volume ratio to match the
units of the differential volume derivation
σ =
A
V
=
2piRov
piR2ov
=
2
Ro
.
Including the reaction for tracking the average island size (4.5), the time-dependent
distributed deposition system is modeled by the following equations
∂
∂t




















[Cu3I3]
[I2]
[Cu]
[CuI]
L
[Cu2]




















= −v
∂
∂ζ




















[Cu3I3]
[I2]
[Cu]
[CuI]
0
[Cu2]




















+




















1 0 0 0 0 0
0 1 0 0 σ/2 0
0 2 −2 −σ 0 0
−3 −2 0 0 −σ 0
0 0 0 0 0 1
0 0 1 0 0 0








































r0
r1
r2
r3
r4
r5




















(4.7)
47
where
r5 = δ
(
rCug + r
Cu
s + r
CuI
g + r
CuI
s
)
from Equation 4.5, where δ  1 is a small parameter that represents the ratio of the
deposition time scale to the species transport time scale, as discussed in Chapter 3,
and it can be estimated as δ = reactor residence time/total deposition time.
4.1.4 Application of the reaction factorization
We apply the reaction factorization developed in Chapters 2 and 3 to the reac-
tion system (4.7) to eliminate redundant dynamic modes and to allow 0,1,2 → 0,
transforming the system into the standard singular perturbation form. The new
reaction coordinates resulting from the reaction factorization are given in Equa-
48
tion 4.8.
y = U · c =




















1 0 0 0 0 0
0 1 0 0 0 0
0 0 0 0 0 1
0 2 −1 0 0 −2
0 0 0 0 1 0
3 2 0 1 0 0




















· c
y =




















[Cu3I3]
[I2]
[Cu]2
2 [I2]− [Cu]− 2 [Cu2]
L
3 [Cu3I3] + [CuI] + 2 [I2]




















=




















z0
z1
z2
x0
x1
w0




















(4.8)
Multiplying the system of ODEs in (4.7) by the transformation matrix U yields the
singular perturbation problem
∂
∂t




















z0
z1
z2
x0
x1
w0




















= −v
∂
∂ζ




















z0
z1
z2
x0
0
w0




















+




















1 0 0 0 0 0
0 1 0 0 σ/2 0
0 0 −2σ 0 0 0
0 0 0 σ σ 0
0 0 0 0 0 1
0 0 0 0 0 0








































1
0g0
1
1
g1
1
2
g2
f3
f4
δf5




















49
where
g0 = K0 [CuI]
3 − [Cu3I3]
g1 = K1 [CuI]
2 − [Cu2] [I2]
g2 = K2 [Cu2]− [Cu]
2
f3 = rCug ([Cu]) + r
Cu
s ([Cu])
f4 = rCuIg ([CuI ]) + r
CuI
s ([CuI ])
f5 = rCug + r
Cu
s + r
CuI
g + r
CuI
s
.
Taking the limit i → 0 for i = 0, 1, 2 and δ → 0 to obtain the DAE system
g0 = 0
g1 = 0
g2 = 0
∂x0
∂t
= −v
∂x0
∂ζ
+ σf3 + σf4
∂x1
∂t
= f5
∂w0
∂t
= −v
∂w0
∂ζ
(4.9)
50
We define the slow time scale τ = δt in which the deposited film property evolves
(see Section 3 in Chapter 3) and let δ → 0 to obtain
g0 = 0
g1 = 0
g2 = 0
dx0
dζ
= −
σ
v
(f3 + f4)
dx1
dτ
= f5 (∗)
dw0
dζ
= 0
(4.10)
The separation of time scales is evident in this formulation from the following two
observations: 1) there is one ODE (labeled *) in time, and 2) the pseudo equilibrium
manifold is defined by a combination of equilibrium relationships and the gas phase
transport/deposition equation. Furthermore, the conserved quantity w0, which is
the conservation of iodine across the reaction in the gas phase, is revealed through
the reaction factorization
dw0
dζ
=
d
dt
(3 [Cu3I3] + [CuI] + 2 [I2]) = 0 ⇒ w0 = wo
where wo = Uwco is the flux for the total number of iodine atoms at any point
within the reactor at any point in time t > 0.
4.1.5 Integration of the DAE system
Before integrating the DAE system (4.10) resulting from the reaction factor-
ization, we must pay special attention to the boundary condition at the start of
51
the deposition zone, well known as one of the main difficulties when solving DAE
systems [27, 29]. Our initial and boundary conditions must not only satisfy the
algebraic equations and the iodine species conservation relationship, but they must
also satisfy boundary condition of pure convective flow at the inlet of the reactor
for the gas phase transport equation in variable x0,
dx0
dζ
∣
∣
∣
∣
ζ=0
= 0 ⇒ 2 [I2]0 − [Cu]0 − 2 [Cu2]0 = 0 ⇒
2 [I2]0
[Cu]0 + 2 [Cu2]0
= 1
Note that this boundary condition in the new reaction coordinates preserves the ratio
of total amount of Cu to I on CuI (the specified inlet condition) at the beginning
of the deposition zone. Assuming the initial condition of an empty reactor, i.e.
c (ζ, 0) = 0, we can now integrate the system (4.10) in time and space in the manner
described below.
First the specified inlet conditions co are transformed into the new coordinate
system to obtain yo
yo = U · co .
Then we project these inlet conditions onto the equilibrium manifold Q, to find the
boundary condition at c0 by solving the following set of equations








Uxc0
g(c0)
Uwc0








=








xo
0
wo








.
where Ux = [2 -1 0 0 0 -2] and Uw = [ 3 2 0 1 0 0] are the rows of U corresponding
to the new reaction coordinates x , and w , respectively. As noted in Chapter two we
recall that even though the manifoldQ is defined in the original reaction coordinates,
52
the above equation yields the projected initial conditions in the full new reaction
coordinates space.
Next we compute the current spatial state of the system by integrating the
new reaction coordinate associated with the slow mode, x0, in space. To do this we
follow an Euler implicit scheme to approximate the difference in space.
dx0
dζ
=
x i+10 − x
i
0
∆ζ
where
x0i = Ux0 · c
i
is the state if new coordinate x0 at the current spatial integration step. Then we
find the composition state of the system at the next step in space by solving the
following system of equations using a Newton-Raphson technique.








(
x i+10 − x
i
0
)
/∆ζ − f
(
U−1yi+1
)
g
(
U−1yi+1
)
wo − wi+10








.
After completing the integration in space, we use the new spatial state of the sys-
tem to evaluet f5 and realize a step integration of the ODE in time (labeled *) in
Equation 4.10
4.1.6 Simulation results and discussion
At the beginning of the simulation we assume there is a number nucSite of
nucleation sites containing a single copper atom each to initialize the island growth
53
model and estimate the surface deposition area at every time integration step. Ta-
ble 4.2 shows the parameters for the nominal operation of the CVD hot-wall reactor
used for experiments in Appendix A and Chapter 6. Figure 4.7 shows the spatial
Parameter Value Units
v 0.03 m/s
K0 4.72×1011 nm6
K1 1.23×10−14
K2 1.09×10−7 nm−3
Ro 0.0195 m
T 750 oC
[CuI ]o 1.48×10−05 nm−3
ηCu 5×10−03
ηCuI 5×10−04
mCu 1.05×10−25 kg/atom
mCuI 3.16×10−25 kg/molecule
nucSites 4.15×10−7 #/nm2
Table 4.2: Parameters used in the simulation of CuCVD at nominal conditions 1 atm at 750oC .
The boat in Figure 4.1 is loaded with 0.1 grams of copper iodide from which [CuI ]o is estimated.
concentration profiles of the precursors in the reactor. It can be seen that the con-
centration of the copper-containing species decreases as the deposition progresses
along the length of the reactor, whereas the concentration of iodine molecules in the
gas phase rises sharply driven by reactions r4 and r1. The concentration axis on the
54
plot in the bottom part of Figure 4.7 is in logarithmic scale to show species with
very low concentrations.
The average size of the growing islands in the film is shown in Figure 4.8.
At the beginning of the ten minute deposition period and in the area closer to the
precursor boat, the concentration of the precursors is higher and thus more copper
is deposited on the surface. The islands consequently grow larger in this zone.
Downstream, the gasp phase concentration of precursors decreases at a reduced
rate and thus the islands that form have a more uniform average size through the
remainder of the reactor.
In Figure 4.9, we compare the direct (top) and surface diffusion (bottom)
deposition rates of copper and copper iodide. Both rates exhibit the same decreasing
trend due to the depletion of precursor species along the reactor. The depleted
precursors result in smaller islands and thus less growth area, which in turn drives
the deposition reactions to slower rates.
55
0 5 10 15 20
rxr position (cm)
1
0
1
2
3
4
5
6
7
8
co
nc
 ×
 1
0 
6
, n
m
-3
Cu3I3
I2
Cu2
CuI
Cu
0.0 0.5 1.0 1.5 2.0 2.5
rxr position (cm)
101
10-5
10-10
10-15
10-20co
nc
 ×
 1
0 
6
, n
m
-3
Cu3I3
I2
Cu2
CuI
Cu
 i i  ( )
co
nc
 ×
 1
0 
6
, n
m
-3
I
I
I
. . . . . .
 i i  ( )
co
nc
 ×
 1
0 
6
, n
m
-3
I
I
I
Figure 4.7: Concentration profiles of precursors along the reactor. The spatial integration of the
system is done on the deposition zone. This zone is the shaded area in the reaction schematic
of Figure 4.1. The boundary condition at the beginning of the deposition zone is indicated by the
filled circle in the copper iodide profile. The light blue shaded area is the physical position of the
substrate on the reactor coordinate.
56
0 5 10 15 20
rxr position (cm)
102
103
L(
nm
)
Lexp  = 1.7 × 10 3
Figure 4.8: Average island size in the growing film along the reactor dimension. The light blue
shaded area is the physical position of the substrate in the reactor coordinate. Note that the model
predicts the average size length of the islands to match the experimental average island size in the
substrate.
57
0 5 10 15 20
rxr position (cm)
102
103
C
u 
at
om
s/
cr
ys
ta
l s
id
e
LexpCu  = 1.7 × 10 
3
0 5 10 15 20
rxr position (cm)
10-2
10-1
100
101
102
103
ra
te
, n
m
-2
s-
1
rg
Cu
rg
CuI
0 5 10 15 20
rxr position (cm)
10-89
10-86
10-83
10-80
10-77
10-74
10-71
10-68
10-65
10-62
10-59
10-56
10-53
ra
te
, n
m
-2
s-
1
rs
Cu
rs
CuI
 iti  
C
u 
at
om
s/
cr
ys
ta
l s
id
e
  .    
 iti  
-
-
ra
te
, n
m
-2
s-
1
I
 iti  
-
-
-
-
-
-
-
-
-
-
-
-
-
ra
te
, n
m
-2
s-
1
I
Figure 4.9: A Snapshoot at time = 5 minutes of the deposition rates of copper and copper iodide
on the surface of growing islands. The plot on the left compares direct deposition rates from the
gas phase, and the plot on the right compares deposition rates by surface diffusion, The scale of
the rate axes is logarithmic.
58
Chapter 5
Numerical tools
This chapter describes numerical tools created to perform the absolute rate
theory calculations used to estimate the adsorption isotherms and surface deposition
rate constants of Chapter 4. The molecular properties of the species and the methods
used to calculate partition functions are structured in object-oriented classes in the
programming language Python and are described in the following sections.
5.1 Chemical species class
A root class called species was created to store the universal constants Avo, kB , h
(Avogrado number, Boltzmann and Planck constants, respectively) to make them
available for inheritance by sub-classes of specific species. While the library of classes
continues to grow, the current classes can divided into two: 1) classes for gas phase
species and 2) classes for adsorbed species.
The classes for the gas phase consist of mono-atomic and polyatomic (including
diatomic) classes. The properties required to instantiate an object of the mono-
atomic gas phase class are simply the molecular weight MW of the gas and the
electronic energy modes. The electronic energy modes are obtained from tables of
electronic energy states that list the energy level ξel and the degeneracy ωel . The
degeneracy is the number of different states (configurations) that have the same
59
energy level. This class contains methods to compute the partition function for the
translational (trans(T)) and electronic (elect(T)) degrees of freedom at a given
temperature T. Details of the derivation and calculation of the partition functions
for this and all classes can be found in statistical thermodynamic references such
as [33, 34].
SPECIES
0i-1+N.1+?
ऊଌഎ༐?
monoatomG
ሓ1+
ᐕถଐ? ??ഘ? ᐐᜑ ?ตᐚ ?ᘙဗ? ᰜᰜᰜᰜᰜᰜᰜ0ဗ1ᴑ?ဗ??ဗ ?‐?1ᴑℐ?1ᴑ
atomSurf
ሓ 
polyatomG
ሓ
 ? ?
1+⌢?ଃᬕᤘ 1+
? ᔎᘋဗ? ᔃᐎᐐ? ?i?␕ฐ? ?ᔃᐎᐐ? ?ఎ?? ᨛᘙ? ᜑ ᰜᰜᰜᰜᰜᰜᰜ0ဗ1ᴑḐᜑἐᜑ‐?1ᴑ?ဗ1ᴑ
adatom
?ล? ሓ?ሓ? ?-?-
?ᔎ?ଐ? ?i?␕?ဗ ?᠍?ᤔ ဗ? చ?iဗ?ఎᔔᨛᘙဗ?
admolec
?ล? ሓ?ሓ? ?-
?-
?ᔎ?ଐ? ?i?␕?ဗ ?᠍?ᤔ ဗ? చ?iဗ?ఎᔔᨛᘙဗ?
[ X^ ]
[ X^ ][ X^ ]
Figure 5.1: Representation of the species classes. atomSurf is a simple class that represents a
surface site, it requires only the molecular weight of the species, its density, and it calculates the
site density [Xˆ ].
60
Gas phase species
The class for polyatomic gas phase species requires the vibrational frequencies
(νi) the symmetry number σF and the moments of inertia Ii of the molecule (or
rotation numbers if available in the CCCBDB/NIST [35] data base), in addition
of the electronic energy levels, to be instantiated. If information to calculate the
rotational energy modes is not available, the species class has another method to
calculate the moments of inertia from the relative location of the atoms to each
other in the molecule [36]. The polyatomic species class has methods to compute the
partition functions for the main energy modes: transitional (trans(T)), rotational
(rota(T)), vibrational (vibra(T)) and electronic (elect(T)).
Surface species
The classes for surface species that represent an adsorbed atom and an ab-
sorbed molecule, are adatom and admolec, respectively. They have common prop-
erties: MWa (molecular weight of the adsorbed species), MWS (molecular weight of
the surface material), ρ density of the surface material and Do the energy of the bond
between the species and the surface. Both species classes have methods for com-
puting the translational, vibrational, rotational and electronic partition functions.
The vibrational mode perpendicular to the surface is estimated within the vibra-
tional methods in both classes. For the admolec class, if the adsorbed molecule has
vibrational modes other than the perpendicular, the frequencies ν’s of these extra
vibrational modes need to be supplied during the instantiation. The partition func-
61
tion calculations of these classes are used to estimate adsorption isotherms through
absolute rate theory as discussed in the next section.
Figure 5.1 shows an schematic describing the different classes and subclasses
for different the types of species in every block rectangle. For a particular block the
first section under the name of the class shows the molecular properties required to
define an object of such class. Then the second block under the name, shows the
functions (methods) that such species supports. e.g. partition function calculations
and thermodynamic state functions such as Helmholtz free energy A(T,P), internal
energy U(T), enthalpy H(T), entropy S(T,P) and G(T,P). This thermodynamic state
functions are calculated with the partition function methods, and where created to
test the validity of the partition function methods in each species in the gas phase.
5.2 Adsorption isotherms from absolute rate theory
To estimate the adsorption isotherm and growth rates we follow Laidler et.
al. [32] who applied the absolute rate theory to develop expressions for adsorption
isotherms. Following their derivation, consider the deposition process illustrated in
Figure 5.2 where a molecule of gas-phase species A deposits on a surface site X and
forms the adsorbed species AX .
The equilibrium criterion at constant temperature yields the following expres-
sion
Keq =
[AX ]
[A] ([X ]− [AX ])
=
NAX/S
(NA/V ) (NX/S)
=
ZAX
(ZA/V )ZX
with units m3 (5.1)
where ZAX , ZA and ZX are the partition functions of the adsorbed species, the
62
+ →
A X AX
Figure 5.2: Schematic of an adsorption reaction. A is the precursor
molecule in the gas phase, X is a surface site available for deposition
and AX is the adsorbed molecule A on the site X .
molecule in the gas phase, and the surface site, respectively. [A] = NA/V is the
concentration of the gas-phase species per unit volume. [AX ] = NAX/S is the con-
centration of surface sites occupied by the adsorbed species and [X ] = NX/S is the
number of surfaces sites available for deposition. S is the area of a surface deposi-
tion site. Note that if [Xˆ ] (sites/area) is the surface site density of the substrates,
S = 1/[Xˆ ].
5.2.1 Copper adsorption isotherm
To make this discussion more concrete, consider the deposition of gas-phase
atomic copper on a copper or quartz substrate. The partition function of a copper
atom in the gas phase will have only ground state electronic, ZA,0 and translational
contributions, ZA,trans . Thus its partition function is
ZA/V = ZA,transZA,0 = (2pimkBTh2 )3/2 exp(−ECukBT
) (5.2)
The adsorbed copper atom can move in only two directions and so loses one
translational degree of freedom relative to the gas phase. However, there is an extra
63
vibrational degree of freedom perpendicular to the substrate. Hence its partition
function is as follows
ZAX = ZAX ,transZAX ,vibZAX ,0
ZAX =
(
2pimkBT
h2
)
·
(
exp (hνo/2kBT )
exp (hνo/kBT )− 1
)
exp
(
−ECu/SiO2
kBT
) (5.3)
where νo is the vibration frequency of the bond between the adsorbed atom and the
surface. According to [32] it can be estimated by
νo =
1
2pi
√
2piro
m∆3
(5.4)
where  is the adsorption energy, ro= 2.35 ×10−10 meters is the bond length for
Cu-Cu [31] and ro = 2.46×10−10 meters for Cu-SiO2 [30], and ∆3 is the volume of
the adsorbed species. This volume was estimated from the surface site density, [Xˆ ],
as ∆3 =
(
1/[Xˆ ]
)3/2
. With these data the vibration frequencies for copper on copper
and quartz substrates are, νo ≈ 2×1012s−1 and ν1 ≈ 5×1011s−1, respectively.
For the surface site, in addition to the electronic ground state, ZX ,0, there is a
vibrational contribution commonly represented with a single frequency ν1, which can
be estimated from the adsorption energy using an equation similar to Equation 5.4 to
be 5.278× 1011s−1 which makes the contribution ZX ,vib = 2.11. Hence the partition
function value of a surface site is given by
ZX = ZX ,vibZX ,0 = 2.11× exp
(
−ESiO2
kBT
)
(5.5)
After applying the above considerations, Equation 5.1 transforms into
[AX ]
[A] ([X ]− [AX ])
=
ZAX ,trans
ZA,transZAX ,vib
ZAX ,0
ZA,0ZX ,0
64
and we have that θAX = [AX ]/[Xˆ ] is the fraction of surface sites occupied with an
atom from the gas phase. Using the ideal gas for [A] = PA/kBT and the equations
for the partition functions we obtain
θAX
1− θAX
=
PA
[Xˆ ]
h
(2pim)1/2 (kBT )
3/2
(
exp (hνo/2kBT )
exp (hνo/kBT )− 1
)
exp
(
−Cu−SiO2
kBT
)
(5.6)
where h and kB are Planck’s and Boltzmann’s constants, respectively. m is the
mass of the atomic species being adsorbed. 1 is the adsorption energy which is the
difference between the energies of the species in the separate and adsorbed state.
Now, if we define the adsorption equilibrium constant
K (T ) =
h
Xˆ (2pim)1/2 (kBT )
3/2
(
exp (hνo/2kBT )
exp (hνo/kBT )− 1
)
exp
(
−1
kBT
)
(5.7)
Equation 5.6 becomes the Langmuir adsorption isotherm
θAX =
K (T ) PA
1 + K (T ) PA
. (5.8)
To calculate the fraction coverage of the surface by copper atoms we use Equation 5.8
and the appropriate adsorption energy 1 for each surface available for deposition,
Cu and SiO2. Then we use Equation 5.7 to calculate the fractional coverage of copper
over a SiO2 surface. The resulting adsorption isotherms are shown in Figure 5.3
From Figure 5.3 we observe that copper, as expected from the adsorption
energy data in Chapter 4, is more likely to deposit on itself than on quartz. These
calculations are used in Chapter 4 to estimate the adsorption isotherm constant,K spads ,
and fractional coverage, θsp, where sp = Cu, CuI , in equations 4.3 and 4.4.
65
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5
pCuI×10
4  [Pa]
10-16
10-14
10-12
10-10
10-8
PCuP
vap
Cu
CuI/Cu
CuI/SiO2
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5
pCu×10
4  [Pa]
10-16
10-14
10-12
10-10
10-8
PCuP
vap
Cu
Cu/Cu
Cu/SiO2
Figure 5.3: Fractional coverage, θ, of copper on deposited copper and on a quartz substrate
as a function of the partial pressure of copper in the gas phase. The red vertical lines
indicate the copper partial pressure at the reactor temperature, PvaCu = 0.133× 10
−4 Pa,
and the copper partial pressure at which our reactor operates, PCu = 4 × 10 −4 Pa.
66
Chapter 6
Experimental deposition and image processing analysis of copper
films
This appendix considers the chemical vapor deposition of copper using copper
iodide as a precursor by the following reaction
CuI(g)→ Cu(s) +
1
2
I2(g)
in the hot wall tubular reactor sketched in Figure 4.1. The copper is deposited on
quartz substrates and optical microscope images are obtained of the grown films.
These images are digitally processed to estimate the size distribution of the copper
clusters that form the film.
6.1 Experimental reactor system
The reactor system and much of the operating procedures for depositing Cu
are the same as those described in Appendix A for the deposition of CuxO films.
An schematic of the reactor is shown in Figure 4.1. Argon is used as a carrier gas
and the copper iodide is fed in powder form in a ceramic boat to evaporate in the
reactor and being carried downstream by the flow of Ar. The temperature settings
for a typical CVD run in this system are in the range 725 - 775 oC . A calibrated
rotameter gauge is used to control the Ar flow rate. While the reactor is heating to
the desired temperature, the Ar carrier gas (fed to the annular inlet section where
67
the CuI boat is located) flow rate is set to a rotameter setting of 60 (570 sccm)
to flush air from the process tube. This Ar flow is held constant through each
experiment and the flow rate of 570 sccm is used for all experiments described in
this appendix. Deposition takes place for 10 min, after which the reactor is allowed
to cool under Ar carrier gas flow for another 10 min before the boat and substrate
are removed. All deposition runs take place under atmospheric (total) pressure.
6.2 Results and analysis
6.2.1 Experimental island size distribution
The deposited films were visually examined under the optical microscope
(Olympus CKX41 with a Photometrics Coolsnap EZ charge-couple device camera).
The resulting image was later processed to have pixels of one color: white for the
substrate and black for the copper molecules deposited.
Figure 6.1: Reproduction of Figure 4.2. The image on the left is a quartz substrate with copper
deposited(red color), the image on the center is a microscope image of an zone in the substrate,
and the image on the right is the black and white version of the microscope image. Details on the
image analysis process are described in this section and also in the caption of 4.2
68
The initial picture from the microscope had a RGB color format and was fil-
tered through a numerical threshold of 127, which is the middle between black (0)
and white (255), to eliminate noise in the image resulting in the black and white ver-
sion shown in on the right of Figure 6.1 (which is a reproduction of Figure 4.2). The
pixels in white are pixels from the substrate (SiO2 which is transparent) while the
black pixels indicate copper particles. The black and white images were processed
using the ndimage tool from the scipy python library, to identify copper islands
(clusters of black pixels) on the substrate.
The image on the right of Figure 6.1 shows that the deposited copper forms
small islands with accumulation a clear tendency for copper crystals to grow along
the film. This suggests that because of the weak bonding of copper with the SiO2
surface, 0.1 eV, the copper deposited diffuses rapidly along the surface until it en-
counters another surface atom or cluster and sticks to it. The scratches on the sub-
strate provide edge and kink sites of higher coordination where the copper adatoms
are more likely to stick as observed in the figure. This concentration on the scratches
and small size island distribution is considered to indicate that the dominant growth
mechanism is the direct deposition of copper atoms on both surfaces. This would
initially seem to contradict the results on the simulation in Chapter 4. However,
this observation suggests that surface migration is important to provide the initial
adsorbed copper atoms (on the kink sites generated by defects on the substrates
surface) where the islands will grow during the deposition. It is then the growth of
the island that is regulated by the direct deposition. In the simulation we assume
that there are initially nucSites (see Table 4.2) on the quartz surface containing a
69
single copper atom where the islands will begin to grow during the deposition, in
accordance with this experimental observations.
For this analysis, the islands are assumed to have a cubic structure since it
is known that copper crystals have a FCC structure. Consequently the number of
pixels of each island can be used to estimate the volume of the cubic island Vp. The
volume of a cubic pixel converted to meters is
Vm3 = Vp × c
3
pix2m = 2.68× 10
−19m3
where cpix2m = 6.45×10−7 m/pixel, is the conversion factor from pixel to micrometers
obtained from the scale at which the camera was used. Finally the number of copper
atoms in the island is estimated by nCu = Vm3/vCu, where vCu is the atomic volume
of copper. The histogram with the island size distribution is show in Figure 6.2.
6.2.2 Average side length analysis
This subsection presents an analysis to validate the the copper crystal size
estimates. For this, the length of the deposition zone in the reactor (see Figure 4.1)
was calculated using the results of the image analysis with a material balance.
Table 6.1 summarizes the analysis of six image samples. Each sample was taken
from a subsection of Figure 6.1, right, in such a manner to avoid scratched areas
that would potentially bias the calculation of island size.
The results in Table 6.1 show that there are in average 1.275×1023 atoms in
every square meter of the substrate, distributed in islands of side length of 1.704
micrometers. In the table, DepLen is the deposition length in the reactor that would
70
0 1 2 3 4 5
Size µm3 1e 17
0
50
100
150
200
250
300
Isl
an
ds
 co
un
t
Island size distribution
Figure 6.2: Overall island size distribution. The size of the islands was estimated taking
the pixel area as the base area of a cubic island.
take to consume all the initial precursor and is estimated as follows
DepLen =
moCuI
MWCuI
ρCu
2piR2o
(6.1)
where moCuI = 0.1 grams is the mass of precursor used in a deposition and Ro = 1.91×
10−2m is the inner radius of the reactor tube. With these data, the average length of
deposition along the reactor is about 19.7 centimeters, which is in good agreement
with visual observations of deposition in the reactor supporting the assumption of
cubic islands.
71
Table 6.1: Analysis of image samples with a scale of 0.645 µm/pixel .
Sample % pixel Cutotal × 1015 ρCu × 1023 DepLen Image size Side Length
number black/Total atoms atoms/m2 cm µm × µm µm
1 15 2.18 1.11 22.4 116×168 1.76
2 18 3.24 1.37 18.1 141×167 1.93
3 17 2.25 1.30 19.1 110×156 1.67
4 14 0.95 1.07 23.2 88×100 1.70
5 19 1.96 1.39 17.9 92×151 1.60
6 19 0.72 1.41 17.6 54×95 1.55
% pixel = number of black pixels / total number of pixels in the image
Cutotal = total number of atoms in the image (total volume/Cu atom volume)
ρCu = density of copper atoms in the image (total number of atoms in the image / image’s area)
[
atoms/m2
]
DepLen = length of deposition in the reactor needed to consume 0.1 grams of copper iodide [cm]
72
Chapter 7
Conclusions and suggestions for future work
A model reduction methodology based on a Gauss-Jordan reaction factoriza-
tion was developed in this thesis. The reaction factorization successfully identified
the dominant dynamics of complex transport-reaction systems, tested the structural
integrity of the kinetic mechanism considered in the system, and provides physical
insight into quantities that are conserved in the system.
To identify the dominant dynamics, the factorization provided a rational ap-
proach to formulating standard singular perturbation representations of complex
transport-reaction systems preserving the dominant dynamics of the original system.
The factorization achieves these results by transforming the original composition
coordinates into new reaction coordinates that are associated with each individual
net-forward reaction rate to the maximum extent possible, i.e. by diagonalizing the
reaction stoichiometry matrix. Even in the case where the diagonalization of the
stoichiometry matrix was not fully realized, the singular perturbation problem was
still valid and accomplished separation of time scales of the system in the outer
solution.
The reaction factorization tested the structural integrity of the kinetic mecha-
nism by recognizing redundant dynamics if they existed in the model; consequently,
avoiding wasteful computations. Furthermore, the factorization revealed conserved
73
quantities in the new reaction coordinates it produced, such as conserved total num-
ber of atoms through the length of the tubular reactor described in Chapter 4.
A major contribution of this thesis was the formulation of the reaction factor-
ization in a form suitable for application to dynamic, distributed-parameter reaction
systems. For these systems, the factorization proved useful not only to separate the
time scales of the system, but as a consequence it also provided a path rigorously
decoupling the time evolution of the system and the spatial distribution so that two
independent ordinary differential equation problems are to be solved rather than
the original set of partial differential equations. Furthermore this approach solved
the problem of how one formulates Danckwerts-type boundary conditions where
gas-phase equilibrium reactions are important.
A series of object-oriented species classes was created to store molecular infor-
mation in each species and structure statistical mechanics computations for abso-
lute rate computations estimations of adsorption isotherms and deposition reaction
rates. These classes have a broad application since they provide an organized way
to automate calculations of fundamental rates that are required in a wide variety of
multi-scale models such as alumina ALD [36].
A crystal growth model was proposed for the growth of copper films on quartz
films at the microscopic level. This model was included in the Cu CVD plug flow
model of Chapter 4. Using this model to track the average size (L) of copper
islands growing on the film. The simulation results predicted the average island size
that matches the experimental average size obtained from image processing of films
Chapter 6 to be in the zone where the substrate is located in the reactor.
74
Finally, the effects of temperature and precursor flow on the composition of
copper oxide films were statistically determined through a design of experiments
and a subsequent surface response analysis of the digitized images of deposited
films. Particularly for the amount of the two oxides, CuO and Cu2O, present on the
films.
Suggestions for future work
1 Factorization of more detailed Cu and CuxO CVD processes.
2 Study of the modeling of 2D materials CVD.
3 Refinement of the film growth model.
Project 1: Factorization of more detailed Cu and CuxO CVD processes
Future work should be extend the reaction factorization process to distributed
CVD, which considers diffusive flux, since these systems are more realistic than the
plug flow model described in Chapter 3. The application is very similar to the
application in Chapter 4, the only difference is that the flux in this case contains
the extra term of diffusion as opposed to only convection in the plug flow model.
75
The original model can be described by the following system of equations
∂
∂t




















[Cu3I3]
[I2]
[Cu]
[CuI]
L
[Cu2]




















=
∂
∂ζ




















JCu3I3
JI2
JCu
JCuI
0
JCu2




















+




















1 0 0 0 0 0
0 1 0 0 σ/2 0
0 2 −2 −σ 0 0
−3 −2 0 0 −σ 0
0 0 0 0 0 1
0 0 1 0 0 0








































r0
r1
r2
r3
r4
r5




















where all the variables have the same meaning as in Chapter 4. Ji = −Di∂[i ]/∂ζ − v [i ] is the flux
of species i , the kinetic expressions are given in Table 4.1, and the parameter for the simulation
are the same from Table 4.2. In matrix form, the system becomes
∂c
∂t
=
∂J
∂ζ
+ S · r.
where S is the reaction stoichiometry matrix and the transformation matrix U that
76
diagonalizes it produces the new reaction coordinates
y = U · c =




















1 0 0 0 0 0
0 1 0 0 0 0
0 0 0 0 0 −2
0 2 −1 0 0 −2
0 0 0 0 1 0
3 2 0 1 0 0




















·




















[Cu3I3]
[I2]
[Cu]
[CuI]
L
[Cu2]




















y = U · c =




















[Cu3I3]
[I2]
[Cu]2
2 [I2]− [Cu]− 2 [Cu2]
L
3 [Cu3I3] + [CuI] + 2 [I2]




















=




















z0
z1
z2
x0
x1
w0




















(7.1)
which allows us to write the system in the standard singular perturbation form
shown below
∂
∂t




















z0
z1
z2
x0
x1
w0




















=
∂
∂ζ




















Jz0
Jz1
Jz2
Jx0
0
Jw0




















+




















1 0 0 0 0 0
0 1 0 0 σ/2 0
0 0 −2σ 0 0 0
0 0 0 σ σ 0
0 0 0 0 0 1
0 0 0 0 0 0









































1
0
g0
1
1
g1
1
2
g2
f3
f4
δf5





















77
Letting 0, 1, 2 → 0 and δ → 0 produces the following semi-explicit DAE system
g0 = 0
g1 = 0
g2 = 0
∂Jx0
∂ζ
= −σf3 − σf4
∂Jx1
∂τ
= δf5
∂Jw0
∂ζ
= 0
(7.2)
which can be solved with a procedure similar to the one described in Chapter 3.
Once the application for the distributed Cu CVD system described above is
completed, the next step should be consider the addition of the oxygen species and
their pertinent reactions with the CuI system. This addition would not change the
methodology for the application of the reaction factorization and can be solved to
study the CuxO system that initialy motivated this research, for which previous
experimental was presented in Appendix A.
Project 2: Study of the modeling of 2D materials CVD.
The island growth mechanism used to estimate the film growth in Chapter
4 can be modified to track the growth of two dimensional islands of different ge-
ometries. This modification enables the factorization in distributed CVD systems
to be applied in the modeling of 2D materials growth that has shown potential for
applications in electronic and optoelectronic devices [37].
78
Project 3: Refinement of the film growth model
Another direction for future work should be the refining of the film growth
model to better track more properties of the film in time and space. For example
a kinetic Monte Carlo (kMC) simulation model enables tracking not only of film
thickness but also of surface roughness [38]. The use of a kMC simulation for the
film growth increases the computational costs of the model. In this context the
model reduction resulting from the reaction factorization can reduce considerably
the computational costs.
An issue commonly encounter in these systems is the noise that the stochastic
results of the kMC introduce on the macroscopic scale transport model. This noise
can be decreased by using a filtering technique, e.g. a Kalman filter, on the kMC
simulation. The resulting reduced model will be well suited to study the potential
of using nanotubes as microreactors such as the ones in [39, 40, 41, 42].
79
Appendix A
Temperature and precursor influence in CuXO film growth and
composition
In this appendix the deposition of Cu2O and CuO by the reactions
CuI + 12O2 → CuO(s) +
1
2I2
2CuI + 12O2 → Cu2O(s) + I2
is considered. A set of experiments was carried out guided by a design of experiments
to study the influence of temperature and precursor flux on the composition of
copper oxide films grown in quartz substrates. The effect of these parameters was
statistically determined with a surface response model analysis.
A.1 Reactor system and substrate preparation
The reactor system consists of a 3.8 cm ID cylindrical quartz process tube
heated by a tube furnace (Fig. A.1). The digital furnace temperature controller
allows time-programmed temperature profiles accurate to 1 oC. The inlet of the
process tube has a custom-manufactured quartz fitting through which the argon
(Ar) carrier and 19% O2 in Ar precursor mixture are separately introduced into
the reactor. Copper(I) iodide is used as the solid-source Cu precursor. Residual
process gas is scrubbed in a chilled methanol bubbler to remove I2 before being
vented through the laboratory hood.
80
Edge si
tअsi
਋
ఈ? ? ??༐ ? ?? e ᐏ?? ? ? ᜖ ?
ᤚ si? e
Figure A.1: CuI CVD reactor process tube and furnace (left); schematic of CVD reactor showing
relative positions of CuI precursor boat, O2 mixture feed tube and substrate location. L1 = 18cm
and L2 = 32cm.
A.1.1 Substrate preparation
The substrates used in this study are 2.54 cm× 2.54 cm square frosted quartz
slides, 2.54 cm wide rectangles of pure copper sheets and ITO-coated glass. Before
each run, the quartz substrates are soaked in a strong HCl solution overnight to
remove previously deposited copper oxide films. Any oxides that remain after the
acid treatment are removed with steel wool and repeated rinsing with acid. Cop-
per substrates are prepared by roughening the surface with extra fine sandpaper;
no preliminary preparation for the ITO-coated substrates was required. All three
types of substrates then are rinsed with distilled water and methanol. Finally, the
substrates are heated to 200 oC in air to dry.
A.1.2 Reactor operating procedure
Temperature settings for a typical CVD run in this system are in the range
of 725 oC to 775 oC. Calibrated rotameter gauges are used to control the Ar and
Ar/O2 mixture flow rates. While the reactor is heating to the desired temperature
the Ar carrier gas (fed to the annular inlet section where the CuI boat is located)
81
flow rate is set to a rotameter setting of 60 (570 sccm) to flush air from the process
tube. This Ar flow is held constant through each experiment and the flow rate of
570 sccm is used for all experiments described in this paper. Before each deposition
run, 0.1 g of the solid copper precursor (CuI) is measured into a ceramic boat; the
boat is cleaned with a strong acid solution after every set of runs. When the reactor
has reached the temperature set point the reactor is loaded first with the substrate,
placed with its leading edge approximately L1 = 18 cm into the reactor from the
entrance collar, followed by the precursor boat containing CuI. The ceramic boat
is placed upstream of the substrate in the zone swept by the Ar carrier gas (see
Fig. A.1). The reactor tube then is sealed with the inlet gas fitting and the oxygen
mixture flow rate is set to the desired rotameter setting. The 19% molar O2 mixture
flow rate corresponding to rotameter settings of M = [5, 10, 15] are [4.9, 6.5, 8.1]
sccm, respectively. Deposition takes place for 10 min, after which the reactor is
allowed to cool under Ar carrier gas flow (no O2 mixture) for another 10 min before
the boat and substrate are removed. All deposition runs take place at atmospheric
(total) pressure.
A.2 Representative films and design of experiments
In this study the effect of varying O2 feed flow rates and reactor temperature on
the deposition mechanism are explored. During the course of this study, the reactor
was operated at three temperatures (725, 750, and 775 oC) and three argon/oxygen
gas mixture flow rates (consisting of rotameter settings of 5, 10, or 15 at each of the
82
three temperatures). This set of nine experiments were repeated multiple times to
collect sufficient data for statistical analysis.
A.2.1 Representative results
Figure A.2 shows representative samples of Cu2O/CuO films produced in our
reactor. Gas flow is from left to right in all cases shown. A striking feature of
these films is the distinct separation between the Cu2O (red) and CuO (black) film
regions. Note that X-ray diffraction (XRD) analysis on samples taken from the each
of the “visually apparent” (red vs. black) cuprous and cupric oxide regions always
shows the presence of a single oxide phase. Furthermore, the films used for this
modeling study were deposited on transparent substrates; visual inspection of each
film revealed identical spatial patterns for the initial and final growth surfaces of
CuO/Cu2O films. Details of the film analysis are reported in Levine[43].
Figure A.2-A shows an increasing fraction of Cu2O from left to right—these
experiments correspond to a decreasing sequence of O2 mixture flow rates. Fig-
ure A.2-B displays a pure Cu2O film deposited on a Cu substrate, produced under
reactor operating conditions that would otherwise deposit a Cu2O/CuO film on a
quartz substrate. Significant differences also were observed in films grown on the
ITO-coated/bare sides of the ITO-coated glass substrates [44]. It is this strong sub-
strate dependency that suggests the predominant role the nucleation and/or surface
reactions play and that motivate the experiments and analysis that follow.
83
A BC
Figure A.2: Cu2O/CuO deposition on quartz (A); O2 mixture flows for the three samples are
M = 60 (left), 30 (center), and 20 (right). Cu2O on Cu substrate is shown as (B) corresponding
to M = 30 O2 mixture flow. The grid shown in the background of A and B corresponds to 1 cm
spacing. Image (C) illustrates the sharpness of the transition between the films; the horizontal
white bar corresponds to 1 mm. This film was deposited with M = 15. For all cases, deposition
was performed at T = 750 oC.
A.2.2 Design of experiments
The effect of varying the reactor temperature T and oxygen mixture flow
meter setting M was investigated. A three-level full factorial experiment design was
created for the two factors T and M ; factor values are listed in Table A.1. The two
factors are scaled by subtracting the nominal operating conditions and dividing the
difference by the maximum deviation used in the three-level design:
m =
M − 10
5
, t =
T − 750 oC
25 K
.
The scaled factors are represented in vector form as x = [m t]T .
Information regarding the nM = 9 sets of experimental conditions were loaded
into a vector of MATLAB paramdata objects. Each object includes measured vari-
able (% cuprous oxide) and the two parameter (factor) names and values; the data
84
M T (oC) m t precursor Cu/O ratio no. exp. Y (% Cu2O) yRSM (reduced) y1D
15 725 1 -1 0.9 3 27.6 27.1 43.7
15 750 1 0 1.2 2 39.5 35.6 64.1
15 775 1 1 1.6 3 31.9 43.9 83.1
10 725 0 -1 1.1 3 37.7 45.9 48.4
10 750 0 0 1.5 9 62.2 54.2 68.6
10 775 0 1 2.0 4 78.5 62.5 87.4
5 725 -1 -1 1.5 3 70.2 64.5 54.4
5 750 -1 0 2.0 5 65.6 72.8 74.4
5 775 -1 1 2.6 3 74.8 81.1 92.9
Table A.1: Experimental conditions defined by a three-level, full factorial design for the two factors
T and M (dimensionless rotameter setting), their normalized values t and m, and the corresponding
measured values for film % Cu2O. The Cu/O ratios correspond to molar values of the atomic
species. The fourth column from the right reports the number of runs performed under each set
of operating conditions; remaining columns give measured Y , empirical yRSM , and 1-dimensional
transport model predictions y1D of film Cu2O%.
structures developed by Leon [45, 46] provide a convenient means of storing the
parameterized data, sorting and manipulating the data, and constructing response
surface models from selected sets of data.
85
A.2.3 Image analysis
Given the large number (35 in total) of experiments that were performed1, a
rapid method for estimating the ratio of Cu2O to CuO over the entire substrate was
required; because of the nature of the deposition reaction and precursor transport
modeling found later in this analysis, direct measurement of film thickness was not
required. Based on the distinct color difference between the oxide phases and owing
to the sharpness of the boundary between them, a computational image-processing
approach was developed to characterize the films.
In our approach, a high-resolution color photograph of each film was read
into a MATLAB image object and was then converted from an RGB to grayscale
image. Histograms of each grayscale image clearly showed a bimodal distribution
in the gray shade distribution (Fig. A.3), with the two peaks corresponding to
the two copper oxide phases. The corresponding pixels of each phase then were
reset to a binary distribution, resulting in the sharply contrasting phases shown in
Fig. A.3; the relative proportions of each phase were determined by simply counting
the number of each of the two pixel values. Experiments corresponding to identical
factor values were averaged, to give the nM = 9 Cu2O percentages shown in Fig. A.3
and in Table A.1.
1Repeated experiments were averaged to give nM = 9.
86
 M =
 15
 
 
 
 
 
M =
 10
 
 
 
 
725o C
M =
 5
750o C
 
775o C
 
0 50 100 1500
1
2
3 x 104 threshold: 70
 
 
m: 1, t: −1
 
 
m: 1, t: 0
 
 
m: 1, t: 1
 
 
m: 0, t: −1
 
 
m: 0, t: 0
 
 
m: 0, t: 1
 
 
m: −1, t: −1
 
 
m: −1, t: 0
 
 
m: −1, t: 1
Figure A.3: The image processing approach to determining film composition from original film
images (top). Note the bimodal distribution in the histogram (center), indicating the fraction of
each form of the copper oxide. Pie charts at the bottom indicate relative ratios of Cu2O (red) and
CuO (black).
87
A.3 Response surface modeling
Given the column vector form of x, we can write the unknown function f that
results in the measured value Y = f (x) +  and a quadratic response surface model
that predicts measurement value as
y = b0 + [b]Tx + xTBx (A.1)
where Y denotes the measured percentage of cuprous oxide in the sample, y is the
corresponding prediction of the cuprous oxide percentage, and  represents measure-
ment error.
In the model (A.1), b0 is a scalar corresponding to the expected value of Y at
the nominal operating conditions, and b2×1 is a vector and B2×2 a matrix of model
coefficients to be determined from the data contained in Table A.1. To see this more
clearly, we can write the model as
y = b0 + b1m + b2t + B1,1m2 + B1,2mt + B2,2t2
From the equation above we can see that element B2,1 is not used and can be set to
zero.
A response surface model of the form (A.1) is identified using the data in
Table A.1, producing a set of coefficients {b0,b,B}. From the discussion above, we
see there are 1 + 2 + 3 = 6 b0, b, and B coefficients to identify in (A.1) from the
nM = 9 data points. Reordering the coefficients into a single column vector c6×1
and using the vector of measurement data Y, the resulting regression problem can
be written as Y = Xc where XnM×6 is a rectangular array constructed by the factor
88
full model reduced model
model coeff vector name value std error sb,k value std error sb,k
b0 c1 61.0 9.93 54.2 3.52
b1 c2 -18.6 5.44 -18.6 4.31
b2 c3 8.30 5.44 8.30 4.31
B1,1 c4 -7.88 9.42 0 0
B1,2 c5 -0.09 6.66 0 0
B2,2 c6 -2.32 9.42 0 0
R2 0.8313 0.7887
σM 13.3 10.5
Table A.2: Model coefficients for the full and reduced response surface models.
values (and all of the relevant products m2, mt, and t2) of the nM experiments.
Having more equations than unknowns, the least squares technique is used to find
the optimal set of model coefficients c = [XTX]−1XTY. This least-squares solution
gives the model coefficients listed in Table A.2.
A.3.1 Analysis of the full model
With our model in hand, we turn to assessing its validity with the main goal
being identification and removal of model terms that are statistically irrelevant. To
begin, using the data contained in Table A.1, we quantify the total variation in the
measured data Y by computing the total sum of squares SST =
∑nM
i=1(Yi − Y¯ )
2 =
3159 where Y¯ is the mean measured value of the modeling data set. The variability
89
not captured by the full model is defined by the sum of squared residuals SSE =
∑nM
i=1(Yi−yi )
2 = 533 where yi is the predicted value of the measured value computed
using the factor values corresponding to measurement Yi . These quantities can be
used immediately to compute the coefficient of determination
R2 = 1−
SSE
SST
= 0.8313
indicating a reasonably good fit to the modeling data. The variance of the modeling
error then is approximated by
σ2M ≈
SSE
nM − p
= 178 and σM = 13.3
where nM = 9 is the number of data points and p = 6 the number of identified
parameters in the model.
These data are used to create Fig. A.4, where the model predictions are com-
pared to the true values; if the model predictions were perfect, all of the blue circles
(modeling data) would line up on the diagonal line. The σM are plotted as red
dashed lines in Fig. A.4. Because of the relatively large value of R2 and that most of
the predictions lie inside ±σM , which allows us to conclude that the model accuracy
is reasonable for the identified response surface model (RSM).
A.3.2 Reduced model
Because of the assumption that parameter estimate errors are normally dis-
tributed with zero mean, we can compute the approximations s2b,k to the true
parameter variance σ2b,k to assess the accuracy of our bi and Bi ,j . These vari-
ances are listed (in terms of standard error) in Table A.2 and are computed as
90
30 40 50 60 70
30
40
50
60
70
measured % Cu2O
pr
ed
ic
te
d 
%
 C
u 2
O
Full model
30 40 50 60 70 80
30
40
50
60
70
80
measured % Cu2O
pr
ed
ic
te
d 
%
 C
u 2
O
Reduced model
Figure A.4: Full model predictions for film composition (left) and the reduced model (right) com-
pared to modeling (blue circles) data; dashed red lines indicate ±σM .
σ2b,k ≈ s
2
b,k = [X
TX]−1k,kσ
2
M where the k , k subscripts denote diagonal elements of
[XTX]−1. The ratio of each estimated parameter standard deviation sb,k and the pa-
rameter value gives some feel for which parameters may be of questionable validity.
This evaluation process is formalized using Student’s t-test whereby the t-values for
each coefficient are first computed as tk = (ck − 0)/sb,k and then used as a test of
whether the coefficient can be statistically distinguished from zero. We then select a
parameter confidence level; in this case we wish to determine which parameters we
are 75% (α = 0.25) sure are actually non-zero. Using the MATLAB function tinv,
we compute tα/2,dof and flag all parameters falling outside −tα/2,dof < tk < tα/2,dof .
The degrees of freedom in this case dof = 9− 6 = 3, and so t0.125,3 = 1.42.
Removing the terms that cannot be distinguished from zero by this test and
recomputing the remaining parameter values using the modeling data set results in
a reduced model predicting the Cu2O percentage y as a function of the two input
91
factors
y = b0 + b1m + b2t (A.2)
with parameter values tabulated in Table A.2. As with the full model, the reduced
model predictions and true measurements are compared in Fig. A.4, where it can
be seen that while there is a bit more spread in the predictions relative to the σM
error bars, the model compares well to the measured data. Furthermore, we find the
coefficient of determination R2 of the reduced model changes only slightly between
the full and reduced models.
Given the relative simplicity of the reduced model, two important trends be-
come clear for this CVD system:
1. Given the negative value of b1 the amount of Cu2O relative to CuO will de-
crease with increasing O2 flow at the nominal operating point, an observation
that makes physical sense based solely on the stoichiometry of the film phases;
2. Because b2 is positive, the amount of Cu2O relative to CuO will increase
with increasing temperature at the nominal operating point, an observation
consistent with Ottosson and Carlsson [9] and with Yoon et al. [47] where
Cu2O/CuO films were deposited by ion beam sputtering in an O2 environment.
Interestingly, this trend is opposite to that reported by Barreca et al. [10] and
Condorelli et al. [48] for MOCVD systems.
92
Appendix B
Derivation of discretization weights
This appendix illustrates the derivation of the discretization weights for the
finite differences of the first and second derivatives used in Chapter 4.
B.1 Second order accuracy Taylor’s expansion
To approximate the first and second derivatives with respect to ζ of a function
f with values defined at three interior collocation points ζi−1, ζi and ζi+1 we use a
the Taylor’s series expansion approximation of f around ζi for the other two points
given by the equations:
fi−1 = fi + f ′i (ζi−1 − ζi ) +
1
2 f
′′
i (ζi−1 − ζi )
2
fi+1 = fi + f ′i (ζi+1 − ζi ) +
1
2 f
′′
i (ζi+1 − ζi )
2
(B.1)
for simplicity in the following derivation we let ∆ζi+1 = ζi+1 − ζi and ∆ζi−1 =
ζi−1 − ζi .
B.1.1 Finite differences of first and second dervatives
To find the discretization weights for the first derivative, we solve the second
equation in (B.1) for f ′′i , and substitute in the first equation to obtain
f ′i =
∆ζ2i+1
ξi
fi−1 +
∆ζ2i−1 −∆ζ
2
i+1
ξi
fi −
∆ζ2i−1
ξi
fi+1 (B.2)
93
where
ξi = ∆ζi−1∆ζ2i+1 −∆ζi+1∆ζ
2
i−1
In a similar manner, to find the discretization weights for the second derivative,
we solve the second equation for f ′′i in (B.1) and substitute in the first equation to
obtain
f ′′i =
∆ζi+1
ξ2i
fi−1 +
∆ζi−1 −∆ζi+1
ξ2i
fi −
∆ζi−1
ξ2i
fi+1 (B.3)
where
ξ2i =
∆ζi+1∆ζ2i−1 −∆ζi−1∆ζ
2
i+1
2
At the lower bound of the interval we cannot use a centered finite difference,
and so we use the values of ζ2 and ζ3 to compute the Taylor series expansion at ζ1
and we find
f ′1 =
(ζ2 − ζ1)
2 − (ζ3 − ζ1)
2
ξ1
f1 +
(ζ3 − ζ1)
2
ξ1
f2 −
(ζ2 − ζ1)
2
ξ1
f3 (B.4)
where
ξ1 = (ζ2 − ζ1) (ζ3 − ζ1)
2 − (ζ3 − ζ1) (ζ2 − ζ1)
2
Finally at the upper bound of the interval (collocation point N), we have
f ′N =
(ζN−1 − ζN)
2 − (ζN−2 − ζN)
2
ξN
fN−2 +
(ζN−2 − ζN)
2
ξN
fN−1 −
(ζ2 − ζ1)
2
ξN
fN (B.5)
where
ξN = (ζN−1 − ζN) (ζN−2 − ζN)
2 − (ζN−2 − ζN) (ζN−1 − ζN)
2
94
Bibliography
[1] Brian G. Willis and David V. Lang. Oxidation mechanism of ionic transport
of copper in Sio2 dielectrics. Thin Solid Films, 467:284–293, 2004.
[2] I. A. Rauf, R. Siemsen, M. Grunwell, R. F. Egerton, and M. Sayer. Gas-phase
nucleation during chemical vapor deposition of copper films and its effect on the
resistivity of deposited films. Journal of Materials Research, 14(11):4345–4350,
1999.
[3] NIki S., Contreras M., Repins I., Powalla M., Kushiya K, Ishizuka S., and Mat-
subara K. Cigs absorbers and processes. Progress in Photovoltaics: Research
and Applications, 18:453–466, 2010.
[4] Maengjun Kim and SungHo Lee SangHo Sohn. Reaction pathways for copper
indium diselenide thin film formation from single multiple InSe/CuSe bilayer
stacks. Thin Solid Films, 531:125–130, 2013.
[5] R. A. J. Shelton. Vapour pressures of the solid copper(i) halides. Transactions
of the Faraday Society, pages 2213–2218, 1961.
[6] R. A. J. Shelton. Vapour-phase deposition of copper by the iodide process.
Transactions of the Faraday Society, 62:222–228, 1966.
[7] W. M. Keeffe. Recent progress in metal halide discharge-lamp research. IEE
PROC., 127:191–189, 1980.
[8] L. S. Hong and H. Komiyama. Chemical vapor deposition of cuox films by cui
and o2: Role of cluster formation on film morphology. J. Am. Ceram. Soc.,
74:1597, 1991.
[9] Mikael Ottosson and Jan-Otto Carlsson. Chemical vapour deposition of Cu2O
and CuO from CuI and O2 or N2O. Surface and Coatings Technology, 78:263–
273, 1996.
[10] D. Barreca, P. Fornasiero, A. Gasparotto, C. Maccato, V. Gombac, and E. Ton-
dello. The Potential of Cu2O and CuO Nanosystems on the Photocatalytic H2
Production. ChemSusChem, 2:230, 2009.
[11] Yatendra S. Chaudhary, Anshul Agrawal, Rohit Shrivastav, Vibha R. Satsangi,
and Sahab Dass. A study on the photoelectrochemical properties of copper
oxide thin films. Internation Journal of Hydrogen Energy., 29:131–134, 2004.
[12] P. E. DeJongh, D. Vanmaekelbergh, and J. J. Kelly. Photoelectrochemistry of
Electrodeposited Cu2O. Journal of the Electrochemical Society, 147:486, 2000.
95
[13] Ryoki Nomura, Yasuharu Seki, and Haruo Matsuda. Preparaton of CuinS2
thin films by single-source MOCVD process using Bu2In(SPr)Cu(S2CNPri2). J.
Mater. Chem., 2(7):765–766, 1992.
[14] D. Arana Chavez, E. Toumayan, F. Lora, C. McCaslin, and R. A. Adomaitis.
Modeling the transport and reaction mechanisms of copper oxide. CVD. Jour-
nal of Chemical Vapor Deposition., 16:336–345, 2010.
[15] U. Maas and S.B. Pope. Simplifying chemical kinetics: Intrinsic low-
dimensional manifolds in composition space, 1992.
[16] S. H. Lam and D. a. Goussis. The CSP method for simplifying kinetics. Inter-
national Journal of Chemical Kinetics, 26(4):461–486, 1994.
[17] M. Hadjinicolaou and D. A. Goussis. Asymptotic solution of stiff PDEs with
the CSP method: The reaction diffusion equation. SIAM J. SCI COMPUT.,
20(3):781–810, 1999.
[18] Aditya Kumar, Panagiotis D. Christofides, and Prodromos Daoutidis. Singu-
lar perturbation modeling of nonlinear processes with nonexplicit time-scale
multiplicity. Chemical Engineering Science, 53(8):1491–1504, 1998.
[19] Nishith Vora and Prodromos Daoutidis. Nonlinear Model Reduction of Chem-
ical Reaction Systems. AIChE Journal, 47(10):2320–2332, 2001.
[20] Marie Nathalie Contou Carrere and Prodromos Daoutidis. Model reduction
and control of multi-scale reaction-convection processes. Chemical Engineering
Science, 63:4012–4025, 2008.
[21] L T Biegler. Differential-Algebraic Equations ( DAEs ). pages 1–40, 2000.
[22] Michael Baldea and Prodromos Daoutidis. Model reduction and control of
reactor-heat exchanger networks. Journal of Process Control, 16(3):265–274,
2006.
[23] S. S. Sujit, A. I. Torres, and P. Daoutidis. Networks with large solvent recycle:
Dynamics, hierarchical control, and a biorefinery application. AIChE Journal,
58:1764–1777, 2011.
[24] Ali Hasan Nayfeh. Introduction to Perturbation Techniques. John Wiley and
Sons, 1981.
[25] R. C. Vieira and E. C. Biscaia. Direct methods for consistent initialization
of DAE systems. Computers and Chemical Engineering, 25(9-10):1299–1311,
2001.
[26] D. E. Schwarz. Consistent initialization for DAEs in Hessenberg form. Numer-
ical Algorithms, 52(4):629–648, 2009.
96
[27] Constantinos C. Pantelides. The consistent initialization of differential-
algebraic systems. SIAM J. Sci. Stat. Comput., 9:213–231, 1988.
[28] U. Mass and S.B. Pope. Simplifying Chemical Kinetic: Intrinsic Low-
Dimensional Manifolds in Composition Space. COMBUSTION AND FLAME,
88:239–264, 1992.
[29] Diana Estevez Schwarz. Consistent initialization for DAEs in Hessenberg form.
Numer Algor, 52:629–648, 2009.
[30] Nuria Lopez, Francesc Illas, and Gianfranco Pacchioni. Ab Initio Theory of
Metal Deposition on SiO2. 1. Cun (n = 1− 5) Clusters on Nonbridging Oxygen
Defects. J. Phys. Chem. B, 103:1712–1718, 1999.
[31] Perla B. Balbuena, Pedro A. Derosa, and Jorge M. Seminario. Dsensity Func-
tional Theory Study of Copper Clusters. J. Phys. Chem. B, 103:2830–2840,
1999.
[32] K. J. Laidler, S. Glasstone, and H. Eyring. Application of the Theory of Ab-
solute Reaction Rates to Heterogeneous Processes I. The Adsorption and Des-
orption of Gases. Journal of Chemical Physics, 8:659–667, 1940.
[33] Stanley I. Sandler. An Introduction to Applied Statistical Thermodynamics.
John Wiley and Sons, Inc, 2011.
[34] Terrell L. Hill. An Introduction to Statistical Thermodynamics. Dover Publi-
cations, Inc., 1986.
[35] NIST Computational Chemistry Comparison and Benchmark Database NIST
Standard Reference Database Number 101 Release 16a, August 2013, Editor:
Russell D. Johnson III http://cccbdb.nist.gov/.
[36] Elizabeth M. Remmers. Master Thesis: Ab initio determination of kinetics for
atomic layer deposition mode. Master’s thesis, University of Maryland, 2015.
[37] Xingli Wang, Yongji Gong, Gang Shi, Wai Leong Chow, Kunttal Keyshar,
Gonglan Ye, Robert Vajtai, Jun Lou, Zheng Liu, Emilie Ringe, Beng Kang
Tay, and Pulickel M. Ajayan. Chemical Vapor Deposition Growth of Crystalline
Monolayer MoSe2. ACS Nano, 8(5):5125–5131, March 2014.
[38] Yiming Luo and Panagiotis D. Christofides. Estimation and control of sur-
face roughness in thin film growth using kinetic monte carlo model. Chemical
Engineering Science, 58(14):3115–3129, 2003.
[39] Hitoshi Ogihara, Sakae Takenaka, Ichiro Yamanaka, Eishi Tanabe, Akira
Genseki, and Kiyoshi Otsuka. Synthesis of SiO2 nanotubes and their appli-
cation as nanoscale reactors. Chemistry of Materials, 18(4):996–1000, 2006.
97
[40] Sy Lee and Ky Choi. Polymerization of Ethylene over rac-Et (1-indenyl)
2ZrCl2/MAO Catalyst Supported on Pseudo-Inverse Opal Silica Particles. In-
dustrial & Engineering Chemistry Research, 51(29):9742–9749, 2012.
[41] Sang Yool Lee, Sung Kyoung Kim, Thao M. Nguyen, Jin Suk Chung, Sang Bok
Lee, and Kyu Yong Choi. Kinetics of styrene polymerization to syndiotactic
polystyrene over metallocene catalyst on flat surface, silica nanotube reactors
and porous silica particles. Macromolecules, 44(6):1385–1392, 2011.
[42] Joong Jin Han, Sang Yool Lee, Sang Bok Lee, and Kyu Yong Choi. Silica Nan-
otube Reactors for Catalytic Polymerization of Styrene and Olefins. Macro-
molecular Symposia, 289(1):25–32, 2010.
[43] R. Levine and R. A. Adomaitis. ISR Technical Report. Technical report,
Unpublished, 2010.
[44] G. Guglietta, T. Wanga, R. Pati, S. Eherman, and R. A. Adomaitis. Solar
Hydrogen and Nanotechnology IV. In Frank E. Osterloh, editor, Proc. of SPIE,
Bellingham, WA, 740807-1, 2009. SPIE.
[45] M. P. Leon. Full Wafer Map Response Surface Models for Combinatorial Chem-
ical Vapor Deposition Reactor Operations m.s. dissertation. M.S. Dissertation,
University of Maryland, 2008.
[46] M. P. Leon and R. A. Adomaitis. Full-wafer mapping and response surface
modeling ttechnique for thin film deposition processes. J. Crystal Growth,
311(13):3399–3408, 2009.
[47] K. H. Yoon, W. J. Choi, and D. H. Kang. Photoelectrochemical properties of
copper oxide thin films coated on an n-Si substrate. Thin Solid Films, 372(1-
2):250–256, 2000.
[48] Guglielmo G. Condorelli, Graziella Malandrino, and Ignazio L. Fragala. Kinetic
study of MOCVD fabrication of Copper(I) and Copper(II) oxide films. Chemical
Vapor Deposition, 5(1):21–27, 1999.
98
Papers
• Arana-Chavez D., Toumayan E., Lora F., McCaslin C., Adomaitis R.A. Mod-
eling the Transport and Reaction Mechanisms of Copper Oxide CVD. Journal
of Chemical Vapor Deposition Vol 16 2010, 336-345.
Presentations
• Arana-Chavez D. and Adomaitis R. A. Chemical Vapor Deposition of Copper:
Study of Gas Phase Nucleation. AIChE Annual Meeting. San Francisco,
California, November 3-8, 2013.
• Arana-Chavez D. and Adomaitis R. A. Modelling spatial variations of cop-
per oxide thin film deposition. ChBE / CHEM / MSE Research Festival.
University of Maryland at College Park, March 2012.
• Arana-Chavez D. and Adomaitis R. A. Modeling Chemical Vapor Deposition
of Copper Oxide Films for Solar Energy Applications and Hydrogen Produc-
tion. SACNAS National Conference. San Jose, California, October 26-30,
2011.
• Arana-Chavez D. and Adomaitis R. A. Modeling Spatial Variations In Thin-
Film Processing of Copper Oxide Films for Solar Energy Applications and
Hydrogen Production. AIChE Annual Meeting. Minneapolis, Minnesota, Oc-
tober 16-20, 2011.
• Dwivedi, V. D. Arana-Chavez, and R. A. Adomaitis Multiscale Simulation
and Optimization of an Atomic Layer Deposition Process in a Nanoporous
Material, Paper 140d 2009 AIChE Annual Meeting, Nashville, TN.
99