# A Self-Consistent Model for Thermal Oxidation of Silicon at Low Oxide Thickness

Authors: Gerald Gerlach; Karl Maser
Journal: Advances in Condensed Matter Physics
 (2016)

---

## Abstract

Thermal oxidation of silicon belongs to the most decisive steps in microelectronic fabrication because it allows creating electrically insulating areas which enclose electrically conductive devices and device areas, respectively. Deal and Grove developed the first model (DG-model) for the thermal oxidation of silicon describing the oxide thickness versus oxidation time relationship with very good agreement for oxide thicknesses of more than 23 nm. Their approach named as general relationship is the basis of many similar investigations. However, measurement results show that the DG-model does not apply to very thin oxides in the range of a few nm. Additionally, it is inherently not self-consistent. The aim of this paper is to develop a self-consistent model that is based on the continuity equation instead of Fick’s law as the DG-model is. As literature data show, the relationship between silicon oxide thickness and oxidation time is governed—down to oxide thicknesses of just a few nm—by a power-of-time law. Given by the time-independent surface concentration of oxidants at the oxide surface, Fickian diffusion seems to be neglectable for oxidant migration. The oxidant flux has been revealed to be carried by non-Fickian flux processes depending on sites being able to lodge dopants (oxidants), the so-called DOCC-sites, as well as on the dopant jump rate.


---

## Body

## 1. Introduction

Thermal oxidation of silicon is one of the most decisive fabrication steps in microelectronics because it allows creating electrically insulating areas which enclose electrically conductive devices and device areas, respectively [1]. Thermal oxidation enables the formation of thin, almost defect-free, and reliable silicon oxide films. This is the reason for the dominance of silicon in semiconductor industry in comparison to other semiconductive materials with better electrical properties than silicon (e.g., germanium of which the first transistors were made [2]). With increasing miniaturization and ongoing downscaling of transistor dimensions, thinner and thinner oxide films are gaining increasing importance because the thickness of gate oxide layers already went down to just a few atoms [3].Thermal oxidation of silicon is based on the migration of oxidant species through the already formed silicon oxide layer and following interfacial reactions. The resulting silicon/silicon oxide interface strongly influences the device properties, for example, of metal-oxide-semiconductor field-effect transistors (MOSFETs). The thinner the silicon oxide film is the more decisive this influence becomes. The silicon/silicon oxide interface comprises a very thin suboxide layer of less than two monolayers and a small region of about three monolayers of distorted silicon in the substrate below the two-monolayer suboxide [4]. At the interface a large number of dangling bonds (unsatisfied valences on immobilized atoms) occur with electronic states within the oxide energy gap. The defect density (about 1010 cm−2) can be reduced by means of annealing by two or more orders of magnitude [4].Deal and Grove derived a very simple formula describing the relation between oxide thicknessL and oxidation time t [5]. Their theory is based on simplified assumptions like linear dependence of the mobile oxidant concentration on position in the oxide and Fickian diffusion with a constant diffusion coefficient. This leads to a linear-parabolic relationship L
=
f
t, where the linear part applies to very short oxidation times and, hence, very thin oxide thicknesses, whereas the parabolic function L
~
t describes the growth at relatively large oxidation times.The simplicity of the model has led to its widespread application, for example, in textbooks [10, 11], even if the authors do not directly refer to it, for example, [12]. However, the model assumptions like constant diffusion coefficient and linear dopant concentration profile do not satisfy the continuity equation on principle. Here, many phenomena occur due to structural differences at the interface between silicon and oxide, for example, stress, dilatation, and relaxation effects. This leads to remarkable differences between the Deal-Grove- (DG-) model and experimentally determined values, for example, published by Blanc [6]. For instance, SIMS analyses of the O-18 isotope in thin silicon oxide layers show exponential concentration decay of O-18 in the depth profile [13]. Reference [14] concludes that the DG-model “provides a successful framework to explain the kinetics of the oxidation process for a wide range of oxide film thicknesses … However, a characterization of the atomic processes occurring during oxidation is well beyond its scope.”Considering the manifold interactions at the silicon-silicon oxide interface on an atomic level, it seems almost impossible to derive a physical model which allows comprising all details in their complexity. For that reason, this contribution uses another approach. Instead of assuming Fickian diffusion we consider the continuity condition, that is, the mass conservation law of the mobile species in the oxide. This reflects that the nature of oxidants keeps constant despite all interactions and non-Fickian effects. In particular, the approach used here comprises the following aspects contrary to other postulates: (i) a non-Fickian oxidant flux in the oxide layer corresponding to the DOCC-sites concept of dopant migration in solids, (ii) a variable reaction coefficient in the first-order reaction of silicon oxide formation at the silicon/silicon oxide interface, and (iii) a consistent formulation for the dependence of the oxidant concentration on depth in the oxide layer. These assumptions have been developed on the base of experimental findings and theoretical studies from previous publications. The present results essentially differ from some conclusions of other statements. The so-called diffusion constant has here been abandoned, whereas the reaction coefficient of the oxidation kinetics appears to be a function of time.


## 2. The Deal-Grove- (DG-) Model

Deal and Grove describe the oxidation process by the following stages of the inward movement of the transported species of the oxidant (Figure1) [5]:(i)
Transport from the bulk of the oxidizing gas to the outer surface where it reacts or is adsorbed: the steady-state fluxF
1 of the oxidant is proportional to the difference of the equilibrium concentration C
∗ of the oxidant in the oxide and the concentration C
0 of the oxidant at the outer surface of the oxide at any given time:(1)
F
1
=
h
C
∗
-
C
0
,
whereh is a gas-phase transport coefficient.
(ii)
Transport of the oxidant across the oxide film towards the silicon: the corresponding fluxF
2 is assumed to be Fickian: (2)
F
2
=
-
D
d
C
d
x
=
D
C
0
-
C
i
L
,
whereD is the effective diffusion coefficient, d
C
/
d
x the concentration gradient in the oxide, C
i the concentration of the oxidant near the oxide-silicon interface, and L the oxide thickness. The right part of (2) follows from the assumed steady-state oxidation condition with ∂
F
2
/
∂
x
=
0.
(iii)
Oxidant reaction at the oxide-silicon interface to form a new layer of SiO2 as a first-order reaction with the reaction rate k: (3)
F
3
=
k
C
i
.Considering steady-state fluxes(4)
F
=
F
1
=
F
2
=
F
3and eliminating C
i and C
0, it yields (5)
F
=
k
C
∗
1
+
k
/
h
+
k
L
/
D
.The growth rate R
=
d
L
/
d
t of the oxide layer follows directly from the flux F considering N as the number of oxidant molecules which form a unit volume of the oxide layer: (6)
R
=
d
L
d
t
=
F
N
=
1
N
k
C
∗
1
+
k
/
h
+
k
L
/
D
.Assuming an initial thickness x
i
=
L
t
=
0, the final oxide thickness L can be calculated from (6) by integration: (7)
N
∫
x
i
L
1
+
k
h
+
k
x
D
d
x
=
∫
0
t
k
C
∗
d
t
,resulting in (8)
L
2
-
x
i
2
B
+
L
-
x
i
B
/
A
=
t
,where (9)
A
=
2
D
1
k
+
1
h
,
B
=
2
D
C
∗
N
,
(10)
B
A
=
k
h
k
+
h
C
∗
N
.Introducing a time τ corresponding to the initial oxide thickness x
i
(11)
τ
=
x
i
2
+
A
x
i
B
,(8) becomes(12)
L
2
B
+
L
B
/
A
=
t
+
τ
.The oxide thickness yields (13)
L
=
A
2
1
+
t
+
τ
A
2
/
4
B
-
1
≈
B
A
t
+
τ
for short time
,
thin oxide,
B
t
+
τ
for long time
,
thick oxide.
B
/
A and B are called the linear and the parabolic rate constant, respectively.Figure 1
Model of Deal and Grove for the oxidation of silicon (after [5, Figure 3]).However, this equation is not self-consistent. If we consider the continuity equation(14)
∂
C
x
,
t
∂
t
=
-
∂
F
2
x
,
t
∂
x
,then it follows from ∂
F
2
/
∂
x
=
0 directly: (15)
∂
C
x
,
t
∂
t
=
0for all 0
≤
x
≤
L. The concentration C
0 of the oxidant at the outer surface of the oxide became time-independent: (16)
C
0
=
C
x
=
0
≠
C
0
tand—given by (4)—all fluxes were F
=
F
1
=
F
2
=
F
3
≠
F
t. Then, (6) gave(17)
d
L
d
t
=
const.or (18)
L
~
t
,which is not observed at long times and thick oxides. Following (2), the concentration C
i of the oxidant near the oxide-silicon interface should decrease with time due to an increasing oxide thickness L
(
t
). This is contradictory to (3) where the flux F
3
~
C
i should remain constant. This shows that the DG-model is not self-consistent [15]. Grove himself agreed to that in a personal communication with one of the authors (Karl Maser): “As the oxide layer grows it will contain more and more oxidant, therefore in reality ∂
C
/
∂
t is larger than 0. This would correspond to a small curvature in the distribution of oxidant instead of it being a straight line. Because of this approximation, strictly speaking, the theory is not internally consistent. However, neglecting this curvature introduces little practical error into the calculation of the rate of oxidation” [16].Summarizing, it can be concluded that the DG-model is limited to the following postulates [5]:(i)
linear dependence of the mobile oxidant concentration on position in the oxide,
(ii)
a constant diffusion coefficient,
(iii)
a constant reaction coefficient in the mass-action law governing the first-order reaction of silicon oxide formation,
(iv)
steady-state conditions [5].Consequently, that model is not able to describe process-inherent changes of interactions and non-Fickian effects in the dynamic system of thermal silicon oxidation.


## 3. Other Models for Thermal Oxidation

Since Deal and Grove have determinedx
i (from the intercept of the linear straight line with the thickness axis) to be about 23 nm, the applicability of the DG-model is restricted to correspondingly higher oxide thicknesses [17]. The DG-model underestimates oxidation rates for thinner oxides significantly. To overcome this discrepancy numerous models were proposed. These models can be distinguished into two classes [10, 17]:(i)
models modifying the DG-model to achieve a better match with experimental data,
(ii)
models assuming particular, different physical mechanisms from the DG-model.The models use either (iii)
empirical or semiempirical expressions describing the underlying mechanisms or
(iv)
generic expressions like power laws of oxide thickness on time.The latter do not presume any detailed knowledge about interactions of dopant (oxidant) migrations and reaction kinetics. However, because they are model-free, they may give indications of new insights into underlying processes.Many models consider parts of the three steps of the DG-model, in particular the transport of the oxidant across the oxide film towards the silicon and the oxidant reaction at the oxide-silicon interface to form SiO2.In [18] a simple extension is proposed in which the diffusion through the amorphous oxide happens via molecular oxygen and the silicon oxidation through the reaction of a small concentration of atomic oxygen. This leads to the following relationship between oxide thickness L and oxidation time: (19)
t
~
1
2
u
+
1
4
exp
⁡
2
u
-
1 with L
~
sinh
⁡
u.The model in [19] assumes two parallel, competing reactions to occur. Molecular oxygen reacts directly with silicon to form both silicon dioxide and atomic oxygen. This is accompanied by the dissociation of O2 as second reaction. The atomic oxygen thus formed either reacts with silicon or recombines to molecular oxygen. The fact that the processes occur in parallel instead of consecutively was explained by the concave curvature of the Arrhenius plot of the observed rate constants [20].Reference [21] describes the oxidation as a process where two—noninteracting—oxidants instead of only one diffuse through the oxide and react at the oxide-silicon interface. Hence, the total oxidation rate yields(20)
d
L
d
t
=
B
1
2
L
+
A
1
+
B
2
2
L
+
A
2with A
1, B
1, A
2, and B
2 as the constants corresponding to (9). This leads to the following oxide thickness versus oxidation time relationship: (21)
L
2
-
x
i
2
+
C
L
-
x
i
-
G
ln
⁡
2
E
L
+
F
2
E
x
i
+
F
=
E
twith(22)
C
=
A
1
B
1
+
A
2
B
2
B
1
B
2
,
E
=
B
1
+
B
2
,
F
=
A
1
B
2
+
A
2
B
1
,
G
=
B
1
B
2
A
1
-
A
2
2
2
B
1
+
B
2
2
.Here, the linear and the logarithmic term of L in (21) describe together the rapid initial oxidation regime. The rate constants B
1, B
2, and B
2
/
A
2 show Arrhenius behavior proportional to exp
⁡
-
E
A
/
k
T with E
A the respective activation energy.Reference [22] (see also [17]) follows the basic assumptions of the DG-model by considering diffusion of molecular oxygen through the oxide followed by a single-step reaction. However, the diffusive flux F
2 of the oxidant is related to the difference in the concentration of oxygen between one interface and another newly forming interface instead of the difference between the surface and the interface oxide-silicon: (23)
F
2
=
F
3
=
k
C
∗
exp
⁡
-
k
D
x
.Finally, the growth law yields (24)
L
=
A
2
ln
⁡
2
B
A
2
t
+
τ
+
1with(25)
τ
=
A
2
2
B
exp
⁡
2
x
i
A
-
1
.The model in [23] (see also [17]) leaves the steady-state assumption for the diffusion of molecular oxygen through the oxide making it consistent with the mass balance at the oxide-silicon interface. As a result of the mass conservation of each species of the system, the oxygen concentration at the oxide-silicon interface becomes almost zero. Solving Fick’s law with the boundary condition that the concentration at the oxide-silicon interface is zero, the growth law yields (26)
d
L
d
t
=
Λ
4
D
1
2
twith(27)
Λ
=
c
O
2
eq
π
c
Si
O
2
exp
⁡
-
Λ
2
erfc
⁡
Λ
.The symbols c
O
2
eq and c
Si
O
2 mean the molar concentrations of O2 in SiO2 at the interface between gas oxygen and silicon dioxide and of SiO2, respectively. Unfortunately, (26) does not lead to similar simple analytical expressions as in the models above.Several authors, for example, the ones of [24], have attributed the enhanced diffusion in very thin silicon oxides to space charge effects. Wolters and Zegers-van Duynhoven [24] compare thermal oxidation of silicon with the classical metal oxidation theory and describe the oxidation process as an electrochemical cell. Oxygen ions move into the silica towards the silicon enabled by thermal energy and driven by the concentration gradient. They charge the silicon negatively and built up an electrical field. Several mechanisms affect the equilibrium of electronic and ionic transport, thus influencing the transport of oxygen and, hence, the oxidation rate. In the case of coupled ionic and electronic currents the kinetics of oxide growth follows a power-of-time law [24–26]: (28)
L
~
t
b
.The exponent b amounts to values between 0.25 (purely parabolic growth) and 1.0.Other authors relate the differences of the oxide growth in thin thermal oxides to compressively strained oxide layers near the silicon oxide-silicon interface which suppress diffusivity [27]. The thickness L
i of the interfacial strained layer is estimated as 0.55
–
1.03 nm for oxidation temperatures between 800 and 1200°C. The inverse growth rate is given as (29)
d
d
L
d
t
d
L
=
N
1
D
C
∗
-
N
1
2
D
2
C
∗
A
+
2
L
d
D
d
L
,where the diffusion constant (30)
D
x
=
D
0
,
0
<
x
<
L
-
L
i
,
D
0
exp
⁡
-
Δ
E
k
B
T
x
-
L
-
L
i
L
i
2
,
L
-
L
i
<
x
<
Lis a function of the depth x. D
0 itself is proportional to exp
⁡
-
E
/
k
B
T. Δ
E is the incremental barrier of the diffusivity at the SiO2-Si interface and k
B the Boltzmann constant. The reaction at the oxide-silicon interface causes a considerably large change in molar volume by a factor of ca. 2.2 which is totally relieved at high temperatures above 900°C and only partially relieved at lower temperatures [28] so that the initial rate constant is deformation-dependent and the parabolic rate constant is stress-dependent.Irene [28] as well as Revesz and Hughes [29] proposed also the transport in micropores within the oxide layer as a reason for the differences between the DG-model and experimental results. Parallel to diffusive flux according to (2), a micropore transport flux F
pores occurs which contributes to the total reaction flux F
3
=
F
2
+
F
pores.Summarizing, all the proposed mechanisms and extensions to the DG-model described above have been strongly doubted by many authors [17, 30]. Up to now, none of these models is widely accepted [10]. In practice, only empirical or semiempirical models are used which introduce additional terms to the DG-model [31]. In [30] it is emphasized that most authors have tried to maintain the validity of the DG-model and have underestimated its limitations. As we have already mentioned above the DG-model is not self-consistent which is agreed on by [24].


## 4. Einstein’s Deduction of Fick’s Law

Fick’s law [32] and the constancy of diffusivities were first theoretically founded by Einstein [33] with the help of a model decisively involving symmetric probability functions of dopant jump distances. Denoting C
x
,
t as the dopant concentration at position x and time t, then its amount will be altered to C
x
,
t
+
τ by jumps of dopant particles which have reached the reference position x from their starting points x
+
X during the time interval t
,
t
+
τ. Here X means the jump distance of the dopant particles. Provided that their probability distribution P is only a function of X but is independent of position, concentration, and so forth, it leads to [34] (31)
∂
C
x
,
t
∂
t
=
x
2
2
τ
∂
2
C
x
,
t
∂
x
2 with(32)
x
2
=
∫
-
∞
+
∞
X
2
P
X
d
X
.Following Einstein, the diffusion coefficient D can be interpreted as(33)
D
=
lim
τ
→
0
⁡
x
2
2
τ
.The constancy of diffusivities is inevitably connected with symmetric probability functions of dopant jump distances, for which Einstein excluded any interactions of these particles with each other. If this assumption does not hold, (31) becomes (34)
∂
C
x
,
t
∂
t
=
x
τ
∂
C
x
,
t
∂
x
+
x
2
2
τ
∂
2
C
x
,
t
∂
x
2
.Since the continuity ∂
C
/
∂
t relates to the divergence of the dopant flux F, ∂
C
/
∂
t
=
-
∂
F
/
∂
x, the first term of the right-hand side of (34) represents a non-Fickian dopant flux which is not only caused by the doping gradient ∂
C
/
∂
x. Accordingly, the quotients x
/
τ and x
2
/
2
τ are variables, which may depend on position, concentration, and so forth. Therefore, dopant transport processes in solids have to be considered in a more generic mechanism than only the Fickian diffusion.


## 5. Non-Fickian Migration of Dopants in Solids: The DOCC-Sites Concept

The DG-model and most of the reported models [18–31] consider only Fickian diffusion regarding (2) as the driving mechanism for the transport of oxygen (molecular, atomic) in the oxide towards the oxide-silicon interface. It is mostly considered a steady-state mechanism, in single cases also non-steady-state [23]. As elaborated in Section 4, the correctness of this assumption has to be doubted, because some dopant fluxes in solids show non-Fickian behavior. For instance, neutron activation analyses of the drive-in process of phosphorus in silicon at 1250°C under pure oxygen ambient atmosphere reveal that [35](i)
uphill-migration of phosphorus in the interior of the silicon crystal takes place near to the Si-SiO2 interface,
(ii)
the total amount of phosphorus contained in the dopant zone of the silicon crystal is independent of the duration of drive-in diffusion. There is not any exodus of phosphorus from the silicon crystal into the SiO2-layer.It has been assumed that the uphill-accumulation is caused by an inhomogeneous distribution of point defects, for example, vacancies, generating a non-Fickian dopant flux component. Non-Fickian phenomena of dopant transport in solids have to be considered on the base of a more generic approach like the concept of occupiable sites for dopant species (dopant-occupiable or DOCC-sites concept, resp.) [34]. This atomistic kinetic framework considers dopant migration in solids as discrete jumps of dopant particles via distinct sites instead of assuming a continuous course of solid-state diffusion. It is based on two conditions: (i)
Sites suitable for occupation by particles of the mobile dopant species must exist—in other words, sites being able to lodge dopants, the so-called DOCC-sites.
(ii)
Particles of the mobile dopant species must have sufficient energy to occupy DOCC-sites.Elementary place change interactions which can occur individually or simultaneously are direct exchanges of atoms and mechanisms involving point defects, for example, vacancy, interstitial, and divacancy mechanisms [36].In [34, 37] the transport of dopant particles via point defects in a cubic lattice is considered in detail. The relationship for the dopant particle flux density comprises four components:(i)
Fickian flux due to the dopant concentration gradientC
′
=
∂
C
/
∂
x,
(ii)
flux driven by the gradientω
′
=
∂
ω
/
∂
x of the dopant jump rate ω,
(iii)
flux driven by the gradientβ
′
=
∂
β
/
∂
x of the DOCC-sites concentration β,
(iv)
drift of particles due to the interior electrical field strengthE:(35)
F
x
,
t
=
α
2
H
·
-
β
ω
C
′
-
β
ω
′
C
+
β
′
ω
C
+
β
ω
C
q
E
k
B
T
.The two contributions in the middle of the parentheses represent non-Fickian parts and are also called crossover components.The symbolα means the lattice constant assuming jumps from one lattice plane to the next. Although silicon oxide is not crystalline in the narrower sense, [38] presents results showing that a transition from the crystalline into the amorphous phase proceeds at the SiO2/Si interface via the crystalline or ordered phase of SiO2. H is the dopant DOCC-site neighborhood constant, which is designated as K in original papers [34, 37]. The dopant jump rate ω describes the probability of the jumping process for a preferential dopant particle to jump into the neighboring lattice plane within a considered time interval. q is the charge of the dopants, T the temperature, and k
B the Boltzmann constant.From this generic model, numerous other effects can be derived like the ones of [39–41] as well as familiar Fick’s law that would follow as a special case for timely and spatially constant parameters β
=
β
0 and ω
=
ω
0. The diffusion constant would yield here D
0
=
H
α
2
β
0
ω
0.In [28] it is stated that Fickian transport of oxygen is not dominant at oxidation temperatures below 1000°C. In [42] it is shown byab initio-modelling of the silicon oxide/silicon interface that thermodynamically stable configurations exhibit defects in the form of three-coordinated silicon atoms, five-coordinated silicon atoms, threefold-coordinated oxygen atoms, and displaced oxygen atoms. Unfortunately, the corresponding parameters of the migration flux parts in (34) are not known in detail for thermal oxidation.


## 6. A Consistent Model for Thermal Oxidation of Silicon

It is the goal of the following considerations to give a summary formulation for the effective oxidant flux (i.e., the fluxF
2 or F, resp.) in the silicon oxide at the interface to silicon. Details of rate process mechanisms, of the oxidation kinetics, and of the nature of the oxidant species will not be considered. We regard the total oxidant flux at the silicon oxide/silicon interface as the decisive parameter determining the growth rate d
L
t
/
d
t of the oxide thickness L
t.Thermal oxidation of silicon is a dynamic process. The stress-strain differences within the oxide as well as between oxide and silicon, which decisively depend on temperature and oxide thickness, alter oxidant migration and reaction kinetics. In this way, the growing oxide is both a reason for and a result of the continuously changing conditions.To develop a model for the thermal oxidation of silicon that is self-consistent, the following assumptions should be made:(i)
non-Fickian oxidant flux contributions in the oxide layer corresponding to the DOCC-sites concept of dopant migration in solids (see Section5),
(ii)
a variable reaction coefficientM in the first-order reaction of silicon oxide formation at the silicon/silicon oxide interface,
(iii)
a consistent approach for the dependence of the oxidant concentration on depth in the oxide layer.(a) Continuity Equation. Due to the assumptions made above, oxidant migration has to be described in a different way than by Fick’s law and the DG-model. For that reason, continuity equation (14) is used to allow the consideration of non-Fickian fluxes of oxidants in the growing oxide.(b) Growth Law for Silicon Oxide. As Blanc outlines in [6] and the references related to Section 2 of this paper show, there have been performed many studies on the rate of growth of oxide thickness. We will here refer to the raw data in [6] for the same reasons as Blanc did. These values are based on oxidation experiments by Hopper et al. [7] who took care to remove the native oxide prior to oxidation. The oxide thicknesses were measured by means of ellipsometry. The measurements cover different oxygen pressure values, temperatures between 780°C and 930°C, and oxidation times starting at a few seconds (corresponding to ca. 1 nm) (Table 1).Table 1
Experimental data from [7] in [6].
Table (series) in [6]
Si orientation
Temperature (in °C)
Oxygen partial pressure (in atm)
Native oxidex
0 (in nm)
Oxidation timet range (in s)
Oxide thicknessL (in nm)
A1
(
100
)
870
1.0
0
18–31,476
0.68–61.8
A2
(
100
)
870
1.61
0
34–18,036
1.37–53.76
A3
(
111
)
780
1.0
0
271–47,838
1.37–36.25
A4
(
111
)
870
1.0
1.2
17–8,310
1.40–38.8
A5
(
111
)
930
1.0
0
30–16,992
3.32–107.3
A6
(
111
)
930
0.57
0
6–11,088
0.87–61.5
A7
(
111
)
870
1.56
0
25–22,569
1.03–96.55
A8
(
111
)
870
1.0
0
35–18,013
1.65–62.23
A9
(
111
)
870
0.56
0
28–24,198
1.43–53.00
A10
(
111
)
870
0.214
0
35–22,019
1.09–30.00
A11
(
111
)
832
1.0
0
58–12,672
1.39–30.05The double logarithmic plots of oxide thicknessL versus time t reveal power law relationships (Figure 2, taken from Figure 1 in [6]) (36)
L
t
=
K
t
b
,where K and b are constants within distinct frame conditions (Table 2). For the example of (111)-Si at 870°C and 1.56 atm, the exponent b amounts to 0.688 (Figure 2). The linearity of the relationship ln
⁡
L
/
nm versus ln
⁡
t
/
s holds down to thickness values of about 5 nm, hence defining the limit of its applicability. Lower thicknesses are difficult to describe not only by the given validity but also due to surface nucleation effects determining the initial stage of silicon surface oxidation up to the first microns.Table 2
Values ofK and b in (36), calculated from the data in [6]. L
m
i
n is the minimum silicon oxide thickness where (36) holds.
Blanc’s series
K in nm/sb
b
Range of validityL
m
i
n (in nm)
A1
0.067
0.645
4.6
A2
0.071
0.659
4.5
A3
0.013
0.744
4.0
A4
0.075
0.691
10.5
A5
0.341
0.569
5.5
A6
0.147
0.649
13.2
A7
0.098
0.688
5.7
A8
0.078
0.682
6.6
A9
0.068
0.659
6.2
A10
0.037
0.669
6.4
A11
0.020
0.782
5.2Figure 2
Oxide thicknessL versus oxidation time t for (a) Si (111) at 870°C and 1.56 atm and (b) Si (100) at 870°C at 1 atm (from [6]).
(a)
(b)The approach of (36) corresponds to similar power-of-time approaches in [43–45]. Reisman et al. [43] found values for b which are neither unity nor 0.5, as might be expected from the DG-model, but exhibited a value of ca. 0.8. Their literature survey presents values of b between 0.25 and 0.8 for dry oxidation depending only on oxidation temperature and oxygen partial pressure. Ngau et al. [45] expanded (36) by a time τ corresponding to the growth of any existing oxide on the silicon: (37)
L
=
K
t
+
τ
b
.In [24, 25] a power series similar to (36) was proposed which was related to coupled ionic and electronic currents in growing oxides (see Section 3). The reciprocal relationship (38)
t
=
τ
+
A
B
L
+
1
B
L
2
+
⋯shows close similarities to the DG-model. Interestingly, it is shown in [25] that 1
/
B and A
/
B are coupled linearly and are also interrelated with τ. It is stated there that this is strongly in favor of models which predict oxidation by power laws.(c) Concentration of the Mobile Oxidant Species in the Oxide. In the following we will use a new generic approach for the concentration of the oxidant: (39)
C
x
,
t
=
C
x
λ
t
,where λ
t is parameter changing by time. The often used ansatz C
x
,
t
=
f
x
/
D
t would—contrary to the given assumptions—lead again to Fick’s equation.Several experiments to determine oxidant profiles in thermal oxides were performed by means of secondary-ion mass spectroscopy (SIMS) using O-16 and O-18 isotopes, respectively [13, 21, 46–48]. The investigations relate to double and triple oxidation steps executed in different isotopic atmospheres. It can be seen in all publications that the resulting profiles differ significantly from Fickian diffusion as is assumed in the DG-model.The measurements in [13, Figure 2] show a clear exponential decay of O-18 up to 40 nm in depth. Extending the second oxidation step from 10 min to 35 min discloses a modified but again exponential profile [13, Figure 3]. The experiments in [21, 47, 48] reveal the O-18 species to be predominantly bonded to lattice sites, so that the interstitial concentration is an insignificant fraction of the migration process. A third oxidation step performed in O-16 after the initial O-16/O-18 oxidation steps shows that the O-18 species stayed with their positions and did not migrate into the new oxide region [21]. The investigations in [46] reveal that the oxide growth obeys a linear behavior up to 50 nm in thickness apart from an initial region of more rapid growth. The O-18 silicon oxide is found to grow also at the outer oxide surface (gas-oxide interface) which corresponds to an uphill-migration of O-18.(d) Oxidant Flux in the Oxide. The oxidant flux F
L at x
=
L
t is given as(40)
F
L
=
F
x
=
L
t
=
N
d
L
t
d
t
,where N
=
4.6
·
10
22
O-atoms
/
cm
3 is the stoichiometric concentration of oxygen atoms of the SiO2 network. Equation (36) as model-free approximation function L
t and the corresponding flux F
L are the essential bases for a consistent model of the oxidation process.To solve continuity equation (14) a function C
x
,
t according to (40) has to be found which fits (36) and (40). The crucial ansatz (41)
C
x
,
t
=
C
S
exp
⁡
-
x
λ
tbeing the basis to develop the present model has been found intuitively by trial and error (Karl Maser).  C
S is a time-independent surface concentration. Bongiorno and Pasquarello calculated in [9] the O2 solubility in amorphous SiO2 and estimated a value for C
S of 19.4
·
10
16 cm−3 for a pressure of 1 atm and a temperature of 1078°C. They found very satisfactory agreement with the corresponding experimental value of 5.5
·
10
16 cm−3 in [8] that corresponds to a deviation of 0.1 eV on the energy scale of the energy minima for the O2 molecule in a-SiO2 which is within the uncertainty of their level of theory.The time derivation of (41) yields (42)
∂
C
∂
t
=
C
S
exp
⁡
-
x
λ
t
·
x
λ
t
2
·
∂
λ
t
∂
t
.Inserting (42) in continuity equation (14) leads to(43)
F
x
,
t
=
-
∫
∂
C
∂
t
d
x
=
-
C
S
∂
λ
∂
t
1
λ
∫
exp
⁡
-
x
λ
·
x
λ
d
xand with ∫
z
·
exp
⁡
a
z
d
z
=
exp
⁡
a
z
a
z
-
1
/
a
2 to (44)
F
x
,
t
=
C
S
∂
λ
∂
t
1
+
x
λ
exp
⁡
-
x
λ
.This formula is valid for the entire oxide layer 0
≤
x
≤
L
t and satisfies continuity equation (14) for the oxidant concentration C
x
,
t. The function λ
(
t
) has to be determined from the boundary conditions at x
=
L
(
t
) and the position of the oxide-silicon phase boundary. Here, (39), (35), and (44) yield:(45)
F
L
=
F
x
=
L
t
=
N
d
L
t
d
t
=
N
·
K
·
b
·
t
b
-
1
=
C
S
∂
λ
∂
t
1
+
K
t
b
λ
exp
⁡
-
K
t
b
λ
.It can be seen easily that a relation of the type(46)
λ
t
=
Q
·
t
bis the sought-after solution of (45). Therefore, one gets the implicit relationship between the constants C
S, N, K, and Q: (47)
N
K
C
S
Q
=
1
+
K
Q
exp
⁡
-
K
Q
.Noteworthily, any connections between parameters named in (47) have been unknown up to the present. Because the values K (from Table 2) as well as C
S and N are known, the ratio K
/
Q is calculable (Table 3). It is much smaller than unity.Table 3
RatioK
/
Q calculated from (47) and the data from Table 2.
C
S
K
/
Q
5
.
5
·
10
16 cm−3 [8]
2.392
·
10
-
6
19.4
·
10
16 cm−3 [9]
8.435
·
10
-
6From (41), (44), and (46) one gets(48)
F
x
,
t
=
b
·
C
x
,
t
·
Q
t
b
-
1
+
x
t
,
(49)
C
x
,
t
=
C
S
exp
⁡
-
x
Q
t
b
.In the previous representation, we intentionally avoided introducing Fickian diffusion as the mechanism for dopant migration as the DG-model does. For that reason, we have now to consider what mechanisms govern the migration processes at thermal oxidation. Equation (49) will be the starting point for these considerations. Inserting x
=
L
t as well as (36) into (49) gives(50)
C
L
=
C
x
=
L
t
=
C
S
exp
⁡
-
L
Q
t
b
=
C
S
exp
⁡
-
K
t
b
Q
t
b
=
C
S
exp
⁡
-
K
Q
,where C
L means the oxidant concentration at the silicon-silicon oxide interface. Putting measured and calculated values of the ratio K
/
Q of Table 3 in (50) gives the decisive result: (51)
C
L
≈
C
S
.In other words, the mobile oxidant concentration throughout the oxide layer has here been revealed to obey the relations (52)
C
x
,
t
=
C
S
=
const
.
,or(53)
∂
C
x
,
t
∂
x
=
0at 0
≤
x
≤
L
t
. This means that the Fickian diffusion is neglectable, so that the flux F
L (45) has to be formulated in terms of non-Fickian components explained in the DOCC-sites concept (Section 5, (35)).Inserting (52) into (48) yields the relation(54)
F
x
,
t
=
b
·
C
S
·
Q
t
b
-
1
+
x
t
,which describes empirically the experimental results. More detailed knowledge about the process kinetics has been deduced from several theoretical investigations [9, 49–52], which disclose neutral oxygen molecules as that species carrying mainly the oxidant flux through the oxide network. Hence, the drift component can be neglected, (55)
q
E
=
0
.The postulate is inferred from numerical analyses, which take into account the difference between the barrier height of atomic oxygen migration and that of molecules. Results of Hoshino et al. [52] based on density functional calculations favor firmly the molecular oxidant migration. The postulate is numerically verified for both α-quartz and β-cristobalite to simulate approximately the silicon dioxide of thermal silicon oxidation. Noteworthily, these studies are more precise than those of Bongiorno and Pasquarello [49] using first-principle calculations and other simulations.The vanishing Fickian oxidant flux following from (51)–(53) gives rise to consideration of the oxidant migration in view of the abovementioned DOCC-sites concept, which takes into account non-Fickian dopant flux components in solids—here oxidant in oxide.Equations (52), (53), and (55) simplify (35) to the relation (56)
F
x
,
t
=
α
2
·
H
·
C
S
·
β
′
ω
-
β
ω
′
,so that the total oxidant flux is carried only by the over-cross components β
ω
′ and β
′
ω. From this together with (14) it follows immediately(57)
∂
F
x
,
t
∂
x
=
-
∂
C
x
,
t
∂
t
=
α
2
·
H
·
C
S
·
β
′
′
ω
-
β
ω
′
′
.Due to (52) and hence ∂
C
(
x
,
t
)
/
∂
t
=
0, (58) gives finally the decisive relation: (58)
1
β
x
,
t
·
∂
2
β
x
,
t
∂
x
2
=
1
ω
x
,
t
·
∂
2
ω
x
,
t
∂
x
2
,which governs the thermal oxidation of silicon in the range of the power law of the oxide thickness on time with respect to (36). It is worth noting that methods to determine β and ω both separately or as a product are hitherto unknown.Assuming a first-order reaction of mobile O2 molecules at x
=
L
t the oxidant flux becomes(59)
F
L
=
F
x
=
L
t
=
M
·
C
x
=
L
t
,where M is the mass-action law coefficient. Inserting again (36), (41), and (46) it yields(60)
F
L
=
M
·
C
x
=
L
t
=
M
·
C
S
·
exp
⁡
-
L
Q
t
b
=
M
·
C
S
·
exp
⁡
-
K
t
b
Q
t
b
=
M
·
C
S
·
exp
⁡
-
K
Q
.Solving (47) for C
S, formulating F
L in terms of (36) and (40), and connecting it to (60) finally lead to(61)
M
=
K
+
Q
·
b
·
t
b
-
1
.Equation (59) can be simplified with the help of (47) being the fundamental relation for the interfacial oxidation reaction in present investigations. Because C
S is much smaller than N, so that the ratio K
/
Q tends to zero, (47) gives the approximate relation(62)
Q
=
K
·
N
C
S
.Therein, C
S, K, and N are known constants. Equation (62) and the inequality K
≪
Q change (61) into the simple relation: (63)
M
=
K
·
N
·
b
·
t
b
-
1
C
S
.
M as the first-order mass-action law coefficient is here a time-dependent parameter which is totally different from the assumptions of (3) for the DG-model.There are hints in the literature that the non-Fickian migration might be caused by stress effects, in particular by nonlinear viscoelastic properties of the silicon oxide which are dependent on the oxide thickness [42, 53, 54]. However, as the results with governing compressive and tensile stresses show, the detailed mechanisms seem still not fully clear [55, 56].


## 7. Conclusions and Summary

In [5] Deal and Grove developed a model for the thermal oxidation of silicon which results in a very simple formula (13). It describes the oxide thickness versus oxidation time relationship with very good agreement for oxide thicknesses of more than 23 nm. Measurements results show that the DG-model does not apply to very thin oxides in the range of a few nm. For that reason, many authors proposed correction terms for DG-equation (13). However, as this paper shows, the DG-model itself is inherently not self-consistent which is also what Grove has agreed on with one of the authors (Karl Maser) [16].The aim of this paper is to propose a self-consistent model that is based on continuity equation (14), not on Fick’s law as the DG-model is.As Blanc’s data in [6] reveal that the relationship between silicon oxide thickness L
(
t
) and oxidation time t is governed—down to oxide thicknesses of just a few nm (cp. Table 3)—by a power law L
t
=
K
·
t
b, where K and b are constants within distinct frame conditions, this has here been confirmed by many other authors, for example, [24–26, 43, 44]. Such a characteristic contradicts the linear-parabolic approach of the DG-model.Differently to gases, dopant migration in solids is not mainly based on Fickian diffusion. For that reason, according to the DOCC-sites concept [34–37], additional contributions have been considered: non-Fickian migration due to spatial gradients both of dopant-occupiable sites and of the dopant jump rate as well as drift due to a gradient of the electric field strength.Other than for Fickian diffusion—where the concentrationC
(
x
,
t
) of dopants is usually described by an ansatz C
x
,
t
=
f
(
x
/
D
t
)—a new generic approach C
x
,
t
=
f
x
/
λ
(
t
) has here been used. This approach matches properly to dopant flux equation (35) of the DOCC-sites concept. Considering the different contributions to the dopant flux, the following conclusion can be drawn:(i)
Given by the calculated [9] and measured [8] values for the time-independent surface concentration of oxidants at the oxide surface and, hence, very low values for the ratio K
/
Q (Table 3), Fickian diffusion seems to be neglectable for oxidant migration. In other words, the oxidant concentration at the silicon-silicon oxide interface is almost the same as the time-independent surface concentration C
S.
(ii)
Inferred from theoretical results about molecular oxygen migration as carrying oxidant species in silicon oxide [9, 49–52], the electrical field strength E in the oxide layer is also neglectable and, hence, so is the drift of oxidant species due to E.
(iii)
As a result, oxidant migration is carried in particular mainly by non-Fickian migration, which is given by the two depth gradientsβ
′
=
∂
β
/
∂
x and ω
′
=
∂
ω
/
∂
x of the concentration β of sites being able to lodge dopants (named DOCC-sites) and of the jump rate ω expressing the probability per time interval that a preferential dopant corpuscle can move from its start position to the adjacent DOCC-site.
(iv)
No methods are known yet, which separately give the parametersβ and ω governing the thermal oxidation of silicon.
(v)
The inequalityK
/
Q
≪
1 yields a time-variable reaction coefficient M of the first-order reaction of the silicon oxide formation at the silicon/silicon oxide interface. This is contradictory to the DG-model.



---

*Source: 7545632-2016-04-19.xml*