Using VICTORY Process to Model the Effect of Silicon Substrate Orientation on Oxidation Kinetics
Introduction
Different crystal planes of silicon are known to have different oxidation rates, e.g. the silicon plane with Miller indices <111> is oxidized approximately 1.7 times faster then the <100> plane [1][2]. During an oxidation process the geometry of a 3D structure may change significantly and the silicon/oxide interface may pass through various crystal planes with different oxidation rates. Therefore, the orientation and type of the silicon wafer affects the resulting geometry of an oxidized structure on this wafer. This paper explains in details how VICTORY Process can be used to model this anisotropic behaviour.
Physical Model
VICTORY Process offers several modes for the simulation/emulation of an oxidation process, namely
 the analytical oxidation mode
 the empirical oxidation mode
 the full physical oxidation mode
The crystal orientation dependence during oxidation is not supported in all modes, but the full physical oxidation mode comprehensively takes care of this effect. Therefore, the oxide growth rate at any point of the siliconoxide interface depends on the direction of the normal to the interface relative to the crystal’s planes when the full physical oxidation mode is applied. This paper deals with the model used in this full physical mode.
The physical oxidation model makes use of two empirical coefficients:
 the linear rate coefficient
 the parabolic rate coefficient
VICTORY Process converts these parameters into microscopical coefficients for oxygen diffusion and and silicon/silicondioxide reaction which are used in the 3D mathematical model. As a rule of thumb,
 for short oxidation time t << τ , the grown oxide thickness h can be calculated according to
(1)
 and for longer oxidation time t << τ, the grown oxide thickness is approximately
(2)
Here
`(3)
and d_{0} is the thickness of initial oxide layer.
The influence of the crystal orientation is modelled by introducing an orientation factor of the linear rate coefficient such that the resulting linear rate becomes
(4)
where is the linear rate coefficient in a reference direction. The factor depends on the direction of the local normal of the silicon/oxide interface.
In the default configuration of the model used in VICTORY Process, the parabolic rate is assumed to be orientation independent which effectively means that only the reaction at the silicon/oxide interface is affected by the surface orientation and not the oxygen transport.
is stored in its functional form in the open material database of VICTORY Process. The data entry <orifac> in the material database which contains the data for is an orientation dependent function. The function has three different values for three different crystal planes, namely <111>, <110> and <100> whereby <111> is the reference direction. Hence refers to this direction. Listing 1 shows the data entry <orifac> in the material database file for the material silicon.
Listing 1: Orientation dependent function <orifac> which models the direction dependence of the linear oxidation rate coefficient. 
This means, for example, that at a silicon/oxide interface point with a normal vector directing towards the <110> crystal plane, the factor equals to 0.833. For an arbitrary crystal direction , not coinciding with any of the <111>, <110> or <100> crystal directions, the orientation factor is defined by an interpolation procedure. The function
(5)
transforms the crystal direction vector into a scalar. The orientation factor is assumed to be a linear function of this scalar passing through the given points:
For an arbitrary direction l the factor is defined by a linear interpolation between these points.
Note, that the material database defines the crystal orientation factor as a function of a normal vector with respect to the basis of the crystal lattice. For simulation purposes the wafer and mask are located and defined with respect to the simulation domain (box). Therefore, the planes of the silicon wafer must be related to the planes of the simulation domain.
Orientation of Silicon Wafer with Respect to Simulation Domain
The crystal orientation of a structure with respect to the simulation domain is set at the beginning of the input deck in the INIT statement by means of the following parameters:
Table 1 shows the correspondence between the simulation domain axis and the crystal planes for different values of the parameters ORIENTATION with parameter SUB.ROT=0. For example for ORIENTATION=111 with SUB.ROT=0 means that the Y axis of computational domain is directed towards the crystal plane with the Miller indices (1,1, 2) and the X axis is directed towards the crystal plane with indices <1 1 0>.


Table 1: Orientation of the coordinate axis of the simulation domain with respect to the silicon crystal for various values of the parameter ORIENTAION when the parameter SUB.ROT=0. 
The statement:
INIT MATERIAL=”Silicon” BOXMIN=”0, 0, 0” \
BOXMAX=”1, 1, 1” \
RESOLUTION=”0.1, 0.1, 0.1” \
INITHEIGHT=0.3 \
ORIENTATION=110 SUB.ROT=43
implies that the vector <110> of a silicon crystal lattice is directed towards the Z axis of the simulation domain, and that the wafer is turned by 43º clockwise around this axis.
Once the mutual location of the simulation domain and the silicon crystal is given, VICTORY Process can determine the orientation correction factor for any point of the silicon/oxide interface. This is done by calculating the normal to the interface point with respect to the basis of the silicon crystal and by interpolating between the values defined for the planes <111>, <110> and <100> according to Equation (5).
Example: Anisotropic and Isotropic Oxidations of Cylindrical Pillar
To demonstrate the dependence of the oxidation rate on the crystal orientation of the substrate, a 3D circular pillar structure as shown in Figure 1 is oxidized. When initializing the process simulation to create this structure the parameters of the INIT statement ORIENTATION=110 and SUB.ROT=0 are set. The structure is oxidized for 5 minutes in a wet oxidation ambient at 1000C. The oxidation simulation is performed twice with two different <orifac> functions:
 anisotropic function (default material data)
 isotropic function (material data are overridden within a user specific local extension of the material database)
Figure 1: Pillar structure to be oxidized. The XY plane of the simulation domain (flat wafer surface) is aligned with the silicon crystal plane <110>. 
The 3D results after oxidation are shown in Figure 2. One can see the difference between the two cases by comparing XZ and XY planes. In the anisotropic case the 3D pictures are not symmetric. To observe the difference in detail we made a 2D cut in the middle of the pillar, parallel to the XY plane for both cases. In case of anisotropic oxidation, the initial circular crosssection is turned into an ellipsoid shape (Figure 3  red lines), while for the isotropic oxidation the silicon/oxide interface and oxide surface keep the form of a circle within a cutplane parallel to the XYplane (plane wafer surface) (Figure 3 green lines).
Figure 2: Pillar structure after 5 min of wet oxidation at 1000 C with an anisotropic oxidation characteristic(top – a)) and with an isotropic oxidation characteristic(bottom – b)). 
Figure 3: The XY crosssection of the oxidized structure. The crystal plane with a direction aligned with axis X is oxidized more due to anisotropy of oxidation. 
For demonstration purposes let us show how the correction factor is calculated by VICTORY Process in this situation. We will find two different correction factors for silicon/oxide interface points with normals aligned with the coordinate axis X and Y. Their normals are =(1,0,0) and =(0,1,0) within the coordinate system of the simulation domain. The coordinates of and in the crystal basis are :
and hence the function values determined by Equation (5) are
Therefore the linear rate correction factors for those points are
It follows that the crystal plane with a direction aligned with the Y axis is oxidized 1.4 times less then the crystal plane with a normal directed towards the X axis. It explains the elongated shape of the silicon/oxide interface shown on Figure 3.
To compare the anisotropic result with the isotropic one we have set an isotropic <orifac> function by creating a local extension of the material database which overrides data from the global material database. Therefore, two folders, namely
smdb/material
have been created in the folder where the experiment is running. The latter folder contains a file
silicon
with the content shown in Listing 2. The environment variable SILVACO_SMDB is set to point to the the first folder. In the local extension of the material database, the values for the <orifac> function were chosen orientation independent and equal to the default value for the <110> plane. Therefore, the resulting oxidation rate along the X direction of the simulation domain remains the same in the anisotropic and isotropic cases.
Listing 2: Orientation independent function in local extension of the material database. 
Conclusion
In this paper we have demonstrated that VICTORY Process is able to take into account the anisotropic properties of oxidation coefficients in threedimensional structures. The appropriate physical model and mathematical procedures behind this supporting software feature have been explained. We have also shown, how to incorporate user’s data (e.g. to create a uniform oxidation rate) into simulations by overriding the relevant sections of the VICTORY Process material database.
References
 B.E. Deal, Journal of Electrochem. Soc, vol. 125, pp. 576, 1978.
 D.W. Hess and B.E. Deal, Journal of Electrochem. Soc, vol. 124, pp. 735, 1977.