Substrate orientation effects on nanoelectrode lithography: ReaxFF molecular dynamics and experimental study

The crystallographic orientation of the substrate is an essential parameter in the kinetic mechanism for the oxidation process. Hence, the choice of substrate surface orientation is crucial in nanofabrication industries. In the present work, we have studied qualitatively the influence of substrate orientation in nanoelectrode lithography using ReaxFF reactive molecular dynamics simulation. We have investigated the oxidation processes on (100), (110) and (111) orientation surfaces of silicon at different electric field intensities. The simulation results show the thickness of the oxide film and the initial oxygen diffusion rate follow an order of (100) > (110) > (111) at lower electric field intensities. It also confirms that surfaces with higher surface energy are more reactive at lower electric field intensity. Crossovers occurred at a higher electric field intensity (7 V nm−1) under which the thickness of the oxide film yields an order of T(110) > T(100) > T(111). These types of anomalous characteristics have previously been observed for thermal oxidation of silicon surfaces. Experimental results show different orders for the (100) and (111) substrate, while (110) remains the largest for the oxide thickness. A good correlation has been found between the oxide growth and the orientation-dependent parameters where the oxide growth is proportional to the areal density of the surfaces. The oxide growth also follows the relative order of the activation energies, which could be another controlling factor for the oxide growth. Less activation energy of the surface allows more oxide growth and vice versa. However, the differences between simulation and experimental results probably relate to the empirical potential as well as different time and spatial scales of the process.


Introduction
Over the past decades, local oxidation-based nanolithography has been proved as an efficient lithographic tool where researchers used AFM tips [1][2][3][4] to fabricate the nanostructures on different substrates. Due to the slowness of the AFM based oxidation lithography, some researchers modified the local oxidation approaches by using other conductive stamps, which facilitates the large-area fabrication [5][6][7]. The principle of this method is shown in figure 1, which is commonly known as Nanoelectrode lithography (NEL). In the NEL system, the pattern transfer is achieved through an electrochemical reaction. When the conductive stamp touches the substrate, a voltage is applied between them. A water bridge is formed between the protrusive parts of the stamp and the substrate [8]. Then the electrochemical reactions occur, which form an oxide film in the touched areas of the substrate. It is commonly represented with the chemical formula: In the past, some simulation studies were carried out to gain an understanding of the oxidation process [9][10][11][12][13]. In our previous work [14], we have utilized a new reactive force field (ReaxFF) developed by the group of Adri van Duin [15] to investigate the molecular dynamics study of NEL process where the structural and morphological features related to the growth of the oxide film have been assessed. However, the orientation-dependent investigation, which is an essential parameter in the kinetic mechanism for the oxidation process, is still missing.
In the case of thermal oxidation, many experimental studies have revealed that the crystal orientation has a dominant effect on the oxidation process [16][17][18][19]. They showed that the orientation of the substrate surface renders a complex effect on the thermal oxidation. Hence, the choice of the crystallographic orientation of the substrate surface becomes more critical in nanofabrication industries. In the present work, ReaxFF MD simulations and experiments have been employed to investigate the substrate orientation effects in the NEL process. We chose silicon substrates with low Miller indices, namely (100), (110) and (111), which are the dominant substrates and structural materials in Integrated Circuits (IC) and device fabrications. Therefore, they have been an excellent model to illustrate some of the fundamental aspects that are involved in the oxidation process. The next sections of this paper are organized as follows. In section 2, the computational and experimental details are described. Section 3 presents the results and discussions, which include the comparison between the substrates with different orientations in terms of oxidation kinetics, charge distribution, oxide film structure, and growth characteristics under different electric field intensities. Finally, in section 4, conclusions are summarized.

Computational and experimental details
In this article, we have utilized the ReaxFF reactive force field, which has previously been successfully employed to study the thermal oxidation [20] and tribochemical wear mechanism [21] of silicon. In order to consider the charge equilibration, the system utilized the charge equilibration (QEq) method [22,23] at each time step. 'Large-scale Atomic/Molecular Massively Parallel Simulator' (LAMMPS) code [24] was used to perform the simulations using the NVT ensemble with the Nose-Hoover heat bath at a time step of 0.1 femtoseconds. The temperature was kept at room temperature (300 K). Verlet integration with a timestep of 0.1 femtoseconds was used to calculate the atom trajectories. OVITO [25] was used to visualize the system. However, the Coulomb energy term in the ReaxFF force field has been changed, as described in [26], which facilitates the incorporation of the applied electric field properly. This modification was made to consider the electronegativity differences and charge conservation, as suggested by Chen and Martinez [27]. We applied periodic boundary conditions (PBC) in both x and y directions while the shrink-wrapped boundary condition was used in z-direction. Simulation models consist of silicon substrates with 10-Å-thick water layers on top of them (shown in figure 2). These water layers correspond to the adsorbed water due to the relative humidity of 90%, which are calculated by using the formula presented in [28]. A reflecting wall was introduced on top of the system to confine the molecules, and an external electrical field was applied between the reflecting wall and the substrate to mimic the presence of a conductive stamp. The details of these models are summarized in table 1.
The experimental observation of the oxide growth for the ultra-short period is unfeasible. However, the oxidation experiment by using a commercial AFM (Bruker DI 3100, USA) with additional circuits to apply voltage pulses will provide insights on the oxidation kinetics of low index silicon surfaces. We have used doped n-type silicon cantilever coated with Cr-Au. The silicon samples were Si (100) with a resistivity of ρ ∼ 10-15 Ω cm, Si (110) with a resistivity of ρ ∼ 1-30 Ω cm and Si (111) with a resistivity of ρ ∼ 1-20 Ω cm. All the nanofabrication and observations were done in the AFM contact mode at room temperature and the relative humidity of 70%.

Analysis of the oxidation kinetics
We have applied an electric field of 7 V nm −1 (i.e. 7 V bias voltage) between the reflecting wall and the substrate for a   duration of 400 ps at 300 K. After the end of the simulation, some adsorption and dissociation of water on the Si surfaces have been observed which finally formed an oxide film (shown in figure 3). The first penetration of oxygen atom to form Si-O-Si bond occurs at 1.1 picoseconds in Si (100) surface, whereas at 1.4 picoseconds in Si (110) surface and 1 picosecond in Si (111) surface (shown in figure 4). From figure 5, it can be noticed that the adsorption and the dissociation of water for the three surfaces have identical trends. At early stages, the reactions take place very fast, after that the oxidation process becomes slow. According to Graf et al, the ionic charges of the oxides affect the charge diffusions in the water bridge that eventually makes the oxidation process slow [29].
The same characteristics can be explained in figure 6, which depicts the evolution of the thickness of the oxide film for the three silicon substrates throughout the simulation. The oxide thickness is determined as the difference between the position (along Z) of the deepest oxygen atom in the oxide film and the position of the silicon atom at the top. For all three cases, there is a quick increase in oxide growth at the very beginning and reaches more than 2.7 Å within 5 picoseconds. Then the oxide thickness shows a steady increase in a layer by layer fashion.
Once the oxygen atom reaches a certain depth, there is no more forward penetration until this layer is saturated. Hence, oxygen atoms fluctuate at this layer until its concentration in oxygen atoms reaches a certain value, after which oxygen atoms start to move to the next layer. It is also observed   that the (110) substrate yields the highest oxide film thickness, whereas the (111) substrate yields the lowest oxide film thickness. At the end of the simulation, the thickness of the oxide films yields the order T(110) > T(100) > T(111). However, this order can be changed at lower electric field intensities that are discussed in section 3.2.
As two types of adsorption (molecular and dissociative adsorption) are observed, one has a significant contribution to the interaction process. From figure 7, it is observed that molecular adsorption contributes more to the interaction process for the substrate (100) and (110), whereas dissociative adsorption dominates for the substrate (111). As shown in figure 8, the three systems have a quick oxidation at the beginning and then a relatively slow oxidation. At the end of the simulations, they exhibit the order (100) > (111) > (110) in terms of the number of oxygen atoms (ML) diffused in the oxide film.
The partial radial distribution functions (RDFs) are useful parameters to characterize the structure of the oxide films. From the partial RDFs shown in figure 9, the first peaks for the mean Si-O bond length have been observed at a position of 1.575 Å, which is very close to the experimental values of~1.61 Å [30][31][32]. After that, some shapes with nondefined orders have been noticed, which suggests the structure of the oxide films is amorphous for all three silicon substrates. Figure 10 shows the oxygen distribution inside the oxide films where the oxygen atoms are counted every 0.4 Å sized bin. The zero position of z refers to the initial surface position and   The charges of the atoms outside the oxide film follow the same characteristics for all three crystallographic orientations. At the outside of the oxide films, the silicon atoms have the charges of around 0 e. The charges for the oxygen atoms (around −0.70 e) and the hydrogen atoms (around +0.30 e) are found very similar to the charges in bulk water. In the oxide zone, the silicon atoms have higher charges as the oxygen atoms surround them. From the charge distribution of the Si atoms (shown in figure 11), it is observed that some of the top-surface silicon atoms in (100) and (110) planes are oxidized with the charges close to 0.63 e, while the top-surface silicon atoms have charges close to 0.45 e in (111) substrate. The average charge of the oxygen atoms in all three oxide films is found as~0.45 e.

Influence of applied electric field
The substrates present interesting features under different applied electric fields. Figure 12 shows the initial oxygen diffusion rates (number of oxygen atoms diffused per second into the oxide layer due to water dissociation), which are calculated for the first 5 picoseconds period. The (100) substrate has the highest initial diffusion rate for either of these three electric fields, whereas the (111) substrate has the lowest initial diffusion rate at lower electric fields (3 V nm −1 and 5 V nm −1 ). At a higher electric field (7 V nm −1 ), the substrate (111) has a slightly larger initial diffusion rate to the substrate (110). Again, there is a crossover in oxide thickness to T(110) > T(100) at 7 V nm −1 , while the other orientation maintains the same order (shown in figure 13). At lower electric field intensities, the order yields T(100) > T(110) > T(111). The initial oxygen diffusion rate and the thickness of the oxide film increase with the increase of the electric field intensity. These types of behaviors have previously been observed for thermal oxidation of various crystallographic orientations of silicon [17][18][19].
To find out a correlation between the thickness of the oxide film and the substrate parameters, we have calculated the surface energies for three different oriented silicon surfaces with the ReaxFF potential function. The results are compared with other results, which indicate a good agreement (shown in table 2). As can be seen that the surface energy exhibits the same order as the order of the oxide film thickness, confirming that surfaces with higher surface energy are more reactive at lower electric field intensity (3 V nm −1 and 5 V nm −1 ).

Experimental results
Few short oxide lines have been generated on each sample by applying 7 V pulses and 3 µm/s of tip velocity. Figure 14 shows the AFM topographic images and line profiles of the fabricated short oxide lines on different silicon substrates.  The oxide height for (110) plane is found to be the largest (2.28 nm), while the height for (100) plane is found to be the lowest (1.49 nm). The (111) plane exhibits an oxide height of 1.73 nm. The same growth order has also been observed for a lower relative humidity of 44% and an applied bias of 12 V. Therefore, ignoring the differences in the resistivity of the substrates, the thickness of the oxide films yield the order T(110) > T(111) > T(100).

Possible explanation for the order of oxide growth
Previous studies on the orientation effect for thermal oxidation of silicon proposed some mechanisms which explain how the crystal orientations affect the oxide growth [17,18,35]. They established a correlation between the initial order of the growth and the density of atoms on planes parallel to the surface. We have observed in section 3.1 that the water molecules react with the silicon at the interface and the oxygen atoms inserted to the Si-Si bonds, which eventually formed Si-O-Si bonds. The reaction sites are the silicon bonds at the surface, and therefore, the oxidation growth is proportional to the concentration of the reaction sites (areal density). Considering the number of reaction sites, the oxidation growth should follow the same order to the areal density of Si atoms, where the (110) plane has the highest areal density of Si atoms, followed by the (111) plane and then the (100) plane (shown in table 1). The activation energy could be another parameter that can influence the oxidation growth. The substrate with less activation energy has more oxidation growth and vice versa. J. R. Lagenza explained the phenomena and investigated the activation energy for different orientations [35]. It is expected that the parallel bonds react easily while the bonds with an angle to the surface react less readily. This is because of a steric hindrance offered by the silicon atom's position in the vicinity of the bond. Consequently, the reaction between the water molecule and a silicon bond depends on their relative positions. As the parallel bonds require less energy to react, the substrate plane with more parallel bonds has less activation energy. The plane (110) has parallel bonds, while the other two substrates have none. Again, the (100) plane has the bonds with a larger angle to the surface plane than that of the (111) plane. Considering these two observations, the activation energies follow the order E(100) > E(111) > E(110).
As the areal density of atoms and the activation energy depend on the substrate orientation, these two parameters could be the main reasons for the variation of the oxide growth. Based on these two parameters, the oxide growth should follow the order T(110) > T(111) > T(100). The experimental results show that the oxide growth follows this order with a small difference between (100) and (111). From the MD simulation results, we have seen a reverse order for these two substrates, while the (110) substrate remains the largest at 7 V nm −1 electric field. However, the differences between simulation and experimental measured values possibly relate to the empirical potential as well as different time and spatial scales of the process. Indeed, the atomic potential does not describe the surface reconstruction of the Si substrate. This would have changed the reactivity (surface energy) of the considered orientations as well as the order of the oxide growth. Further MD studies at longer simulation time will be required to investigate these varieties of behaviors and to study which factors are responsible for these irregularities. Nevertheless, it will require very long computation times.

Conclusions
The effects of crystallographic orientations of substrates on the NEL oxidation process has been investigated for the (100), (110) and (111) surfaces of single-crystal silicon through ReaxFF MD simulations and experiments. Simulation results revealed that molecular adsorption of H 2 O contributes more to the interaction process for the substrates (100) and (110), whereas the dissociative adsorption dominates for the substrate (111). The RDF analysis suggests the oxide film as an amorphous structure for all three crystallographic orientations. This study further demonstrated that the thickness of the oxide film exhibits the order T(100) > T(110) > T(111) at lower electric field intensities. A crossover occurs at higher electric field intensity (7 V nm −1 ) for the substrates (100) and (110). Besides, it is found that there is also a slight crossover in the initial oxygen diffusion rate to R(111) > R(110) at 7 V nm −1 . The initial order at lower electric field intensities was R(100) > R(110) > R(111). The oxygen diffusion rate and the thickness of the oxide film increase with the increase of the electric field intensities. AFM oxidation experiment shows the different order of the oxide thickness for the (100) and (111) with a small difference. However, molecular dynamics simulations were performed for 400 picoseconds only, while experimental investigations show large-scale features. Therefore, further studies are required with a long simulation run time to elucidate some of the many unknowns concerning the origin of the anomalous features. ledge the use of the EPSRC (EP/K000586/1) funded ARCHIE-WeSt High-Performance Computer at the University of Strathclyde.

Data statement
All data underpinning this publication will be available from the University of Strathclyde Knowledge-Base at https://doi.org/10.15129/f92d7c9f-a52f-4908-afbf-c50330be2edd.