# Vibrational Characteristics of Au-Doped Si Nanowires: A Molecular Dynamics Study with a Modified Embedded Atom Method Potential Developed for Si-H-Au Systems

Authors: Bin Liu; Yun’an Zhang; Xinfeng Wang; Yanjiao Li; Yingjuan Yue
Journal: Journal of Nanomaterials
 (2022)

---

## Abstract

Silicon nanowires (SiNWs) are promising structures as resonators for applications in ultrasensitive force and mass transducers. In the present work, large-scale molecular dynamics simulations were conducted to explore the effect of gold (Au), a commonly used catalyst, on the vibrational properties of [1 0 0] and [1 1 0] SiNWs. The force field was established using the modified embedded atom method (MEAM) formalism. The parameterization focused on the properties of the hydrogen-covered SiNWs doped with Au, involving the formation energy as well as the structural and vibrational characteristics. The vibrational characteristics of a clamped-clamped SiNW were examined through excitation using a sinusoidal velocity field after energy minimization and thermal equilibration. Different doping concentrations and distributions can significantly change the frequency response, quality factorQ, and beat phenomenon. The natural frequency f0 decreased with increasing Au concentration, whereas the effect of the impurity distribution on f0 was negligible. Furthermore, Q was sensitive to the Au concentration and distribution, and the increased concentration led to an overall increment in Q, accompanied by considerably increased scattering. Besides, the beat period of [1 1 0] SiNWs showed a strong positive correlation with the concentration. The vibration along the elementary axes cannot produce the beat phenomenon. Deflection of the elementary axes in a very small range would lead to a higher energy dissipation rate. In addition, owing to the considerable difference between the atomic masses of Au and Si, randomly distributed Au atoms significantly disturb the symmetry of the SiNWs. However, the elementary axes of the [1 1 0] SiNWs remained stable, even in the models with extremely radial or longitudinal impurity segregation. In contrast, we failed to capture the distinguishable elementary axes for most of the doped [1 0 0] samples. The features described above indicate that the impurity Au can effectively tune the resonant frequency of the [1 1 0] SiNWs, meanwhile, causes no large frequency shift that could be potentially caused by the deflection of the elementary coordinate system.


---

## Body

## 1. Introduction

The distinctive properties, including surface effects, quantum effects, and electro-optic properties, facilitate nanowires (NW) to be widely applied in critical domains, such as biosensing, thermoelectrics, and new generation lithium batteries. Size reduction leads to outstanding natural frequencies in the range of GHz to THz and excellent energy efficiency; hence, they are used as key structures in ultrasensitive mass [1] and force [2] transducers.The influence of impurities on the vibrational characteristics of NWs is a topic of research interest. Zhao et al. [3] found that rare metal elements can modulate the sensitivity of ZnO nanowires. Islam et al. [4] systematically studied how the guest element distribution governed by grading functions modifies the mechanical and vibrational properties of functionally graded Cu-Ni NWs. Georgakaki et al. [5] investigated the effects of the Si percentage on the natural frequency f0 and quality factor Q in SixGe1-x and Si/SixGe1-x nanowires. To be strict, the minor functional elements in the above studies are not impurities because they are intentionally introduced. To date, we have not found any other reports on the effects of impurities on the vibrational performance of SiNWs. In the fabrication processes of SiNWs, impurities, such as metal catalysts in the bottom-up procedure and high-energy ions in the top-down procedure, are introduced inevitably. Therefore, this study investigated how the dissolving Au [6], commonly used for catalyzing the unidirectional growth of SiNWs [7], modulates the vibrational properties of SiNWs.Studies on the vibrational characteristics of SiNWs mainly focus onf0 and Q. SiNWs vibrate near f0 when functioning as resonators, and a sufficiently high Q is required to ensure the accuracy of force and mass detection. Experimental, analytical, and numerical studies addressing these aspects have emerged in recent years. Technical parameters, including vibration modes and their corresponding frequencies [8] and Q determined experimentally, are used to evaluate the vibrational performance of SiNWs. In addition, experimental measurements can be used to control the cubic nonlinearity of SiNW resonators [9]. On the analytical front, the mathematical models developed based on continuum theory can calculate the f0 of an NW with surface effects [10] and characterize the nonlinear vibration influenced by surface stress [11]. Numerical simulations can accomplish the abovementioned tasks with appropriate initial excitation [4, 5, 12, 13]. Compared with analytical methods, numerical methods, e.g., molecular dynamics (MD), are advantageous because they offer flexibility for structural modeling for nonideal NWs. In other words, NWs with internal or surface defects are in line with physical reality. Obviously, numerical methods are more economical and can also extract vibrational features without external noise disturbance. Besides, time-dependent properties and kinetic relationships are beyond the scope of the former two methods. However, in numerical simulations, the force field used to describe the behaviors of an NW should be sufficiently reliable.In this study, large-scale MD was used to investigate the vibrational characteristics of Au-doped [1 0 0] and [1 1 0] SiNWs. The influences of the impurity concentration and distribution on the frequency response,Q, and beat phenomenon were systematically investigated. This paper is divided into two sections. The first section presents the optimization of the parameter sets for the Si-H-Au ternary MEAM potential and the reliability validation. The second section explains the simulation results and phenomena in detail.


## 2. Parameterization of the MEAM Potential for the Si-H-Au Ternary System

This work employed the second-nearest MEAM (2NN MEAM), and its mathematical scheme is briefly introduced in an earlier Ref. [14]. The 2NN MEAM has been successfully applied to describe metals, semiconductors, and gas elements, and it possesses excellent scalability for multielement systems [15]. It also provides a good description of the Si surface properties [16], which are especially important for this work. We did not choose the state-of-the-art surrogate model [17, 18] fitted by machine learning; although, it is more accurate and has outstanding transferability. More importantly, it is several orders more computationally expensive and thus is not suitable for dynamic processes in nanoseconds. All empirical models are confined to Newtonian classical mechanics, which cannot accurately predict all properties of a particular material. In particular, when the simulations involve bond breaking and recombination, the credibility of the empirical potentials is questionable. Therefore, the SiNWs were allowed to vibrate only in the elastic range.There are several studies on the MEAM potentials for Si-Au binary systems. Ryu and Cai [6] fitted the parametric set for Si-Au MEAM potential based on the binary phase diagram, which was used to simulate the growth of SiNWs catalyzed by Au. Kuo and Clancy [19] optimized the MEAM potential according to several hypothetical Si-Au compounds, and it was applied to simulate the evolution of the Si-Au interface. Empirical potentials can be optimized for a particular simulation scenario. In this study, the vibrational characteristics of Au-doped SiNWs were highlighted. The theoretical expression of the natural frequency of a clamped-clamped NW is given below:
(1)fn=ωn2πL2EIρA,where the length L and moment of inertia I are determined by the size and cross-sectional shape of the NW. ωn is the eigenvalue of the characteristic equation cosωncoshωn=1. From the simplest perspective, impurities modulate the natural frequency by changing the elastic modulus E and density ρ. ρ is determined by the impurity concentration, whereas E is mainly related to elastic constants and lattice distortion. Following the above guideline, the 2NN MEAM potential was optimized in this work. In addition, MD can only deal with thin NWs (<20 nm). The scale that can be handled is smaller owing to the complexity in the formalism of the MEAM. H atoms often terminate the surface of an ultrathin NW after removing the surface oxide with HF [20]. Therefore, the Si-H-Au ternary system is considered to guarantee that our models are in line with engineering specimens.Several versions of the MEAM model have been developed. This work used the common formalism [21] to calculate the energy as follows:
(2)E=∑iFiρ¯i+12∑j≠iϕRij.The first term is the embedding energy for atomi with background electron density ρ¯i. The second term describes the pair interaction, which is the sum over all neighbors j of i within the cutoff distance. The following equation was chosen to calculate the background electron density ρ¯i.
(3)ρ¯i=ρi021+exp−Γi,where ρi0 is the spherically symmetric partial term, and Γi is the angular dependence. The left mathematical description of the MEAM potential is in line with Ref. [15].Table1 summarizes the MEAM parameters of elements Si, H, and Au. The sets for Si and H were extracted from earlier publications [16, 22], respectively.Table 1
MEAM parameters for Si, H, and Au. The units of cohesive energyEc and bulk modulus B are eV and GPa, respectively. The cutoff distance is 6 Å.
Ele.aEcαBβ(0)β(1)β(2)β(3)At1t2t3CminCmaxdattracdrepulρ0Si [16]4.634.90993.552.5007.500.581.805.25-2.611.422.80001H [22]2.2252.960.74c2.963.003.002.502.500.20-0.100.002.002.800020.5Aub3.936.601804.994.684.00-2.191.325.903.2916.462.002.8000.0660.48aThe reference structures of unary Si, H, and Au are diamond, dimer, and face-centered cubic (FCC), respectively. bDeveloped in this work. cThe distance between H atoms in a dimer instead of the bulk modulus.


### 2.1. MEAM Potential for Au

Several available MEAM potentials for Au [6, 19, 22, 23] were not chosen for our work, as the parameters had to be optimized for obtaining better predictions of the fundamental properties of the Si-Au binary system. The properties referred for optimizing the Au unary MEAM potential include lattice constant, elastic constants, surface energy, and formation energy of vacancies. The MEAM values of these properties are satisfactory, as shown in Table 2. The experimental datum of the (1 0 0) surface energy was considered in the fitting. Therefore, the MEAM prediction showed that the energy of the (1 0 0) surface was much lower than that of the (1 1 0) surface. This is inconsistent with the first principle calculations but agrees with the tight-binding calculations. In addition, the MEAM potential provides higher values of the structural energy differences in face-centered cubic (FCC)/body-centered cubic (BCC) and FCC/hexagonal close-packed (HCP) structures compared to the other calculation methods, suggesting that the current MEAM potential provides a more stable FCC structure.Table 2
MEAM predictions of the lattice constanta0 (Å), elastic constants Cij (GPa), formation energy of vacancy (eV), surface energies (eV), and structural energy differences EFCC/BCC and EFCC/HCP (eV), in comparison with experimental values and higher-level calculations. The energies of vacancy formation and surface are obtained after relaxation.
PropertiesExpt.First-principleTight-bindingMEAMOthersThis worka04.07e4.068e, 4.079c—4.073e4.078C11201.63a217b, 187c—201.5b204.75C12169.67a171b, 168c—169.7b167.20C4445.44a47b, 37c—45.4b, 38d42.02B180.32a174c——179.72Ev0.89±0.04d0.55b0.49d0.91b, 0.7d0.89γ(100)1.54b1.97b1.54d1.08b, 0.95d1.55γ(110)—1.10b1.85d1.05b, 1.00d1.83γ(111)—0.92b1.48d0.90b, 0.70d1.42EFCC/BCC——0.019d0.02d0.18EFCC/HCP——0.0086d0.021d0.09a[24], 0 K values extrapolated from smooth curves, which are determined experimentally. bValues calculated in Ref. [23] or quoted therein. cDensity functional theory (DFT) calculations in this work. dValues calculated in Ref. [19] or quoted therein. eValues calculated in Ref. [25] or quoted therein.Next, we performed density functional theory (DFT) calculations to obtain the necessary reference data for developing binary and ternary MEAM potentials. The exchange correlation functional is approximated by generalized gradient approximation (GGA) revised for solids by Perdew-Burke-Ernzerhof [26], with plane wave expansion up to the cutoff energy of 300 eV determined by convergence tests. The electronic structure as a function of the positions of atomic nuclei was tackeled using the projector augmented wave method. The relaxation was repeated until the total energy reached below 10-5 eV and force below 0.01 eV/Å. The Brillouin zone was sampled with the Monkhost-Pack scheme. For calculating the cohesive energy, lattice constant, and elastic constants, 10×10×10k-points were used in the model with 1×1×1 unit and 6×6×6k-points for point defects in 3×3×3 units. The computations were completed using Vienna Ab initio Simulation Package (VASP) [27–29].Some basic physical properties of Si and Au were calculated and compared with the experimental values (given in parentheses) to verify the reliability of the above simulation settings. The lattice constant of the Si crystal was 5.439 Å (5.431 Å), the cohesive energy 4.49 eV (4.63 eV), and the elastic constants C11, C12, and C44 were 157 GPa (166 GPa), 62 GPa (64 GPa), and 75 GPa (80 GPa), respectively. The calculated and experimental values for Au are listed in Table 2.


### 2.2. MEAM Potentials of Si-Au, Si-H, and H-Au Binary Systems

#### 2.2.1. Si-Au Interaction

For the Si-Au binary system, it is necessary to precisely predict the lattice distortion and changes in the elastic and energy properties of bulk Si when Au is introduced. The Si7Au compound was modeled by replacing one Si atom in the cubic unit with an Au atom. This structure is considered for parameter tuning as it represents the bulk Si doped with a high Au concentration of 12.5%. Specifically, the lattice constant, cohesive energy, and elastic constants of Si7Au and formation energies of the substitutional and interstitial Au in bulk Si were used in the parameterization. Owing to the absence of stable compound structures in the Si/Au phase diagram, b1-NaCl was chosen as the reference structure for the Si-Au binary system. Some other reference structures were considered but failed to yield predictions close to the reference values. As it is a hypothetical structure, the lattice constant and cohesive energy of b1-SiAu were adjusted slightly based on DFT data to obtain overall agreement with the reference data. Table 3 lists the MEAM parameters of the binary Si-Au system.Table 3
MEAM parameters of Si-Au, Si-H, and Au-H binary systems. The units of cohesive energyEc and nearest neighbor distance re are eV and Å, respectively. Refer to [15] for the definition of each parameter.
ParameterSi-AuSi-HH-AuReference structureb1-SiAuDiamond-SiHDimer-HAuEc5.653.492.20re1.851.5851.50α2.1332.7743.990Cmin (A-B-A)2.002.202.00Cmin (A-A-B)2.001.002.00Cmin (A-B-B)2.001.002.00Cmin (B-A-B)2.002.002.00Cmax (A-B-A)3.803.502.80Cmax (A-A-B)2.503.802.80Cmax (A-B-B)2.802.802.80Cmax (B-A-B)2.802.802.80dattrac000.18drepul04.25-0.015The performance of the potential for the Si-Au system was examined. First, the differences in the lattice constant, cohesive energy, and elastic constants between Si and Si7Au were measured. The cohesive energy was calculated as
(4)Ec=NSiESi+NAuEAu−EtotNSi+NAu,where ESi and EAu are the free energies of Si and Au atoms, respectively, NSi and NAu are the numbers of Si and Au atoms in the compound, respectively, and Etot is the total energy.The results obtained from DFT and MEAM are compared in Table4. When the concentration of Au is as high as 12.5%, the lattice distortion and energy shift can be accurately determined. The lattice constant of Si7Au increased because of the repulsion between Si and Au atoms.Table 4
Comparisons of the lattice constanta0 (Å), cohesive energy Ec (eV), and elastic constants Cij (GPa) between Si7Au and Si.
a0EcC11C12C44BDFTMEAMDFTMEAMDFTMEAMDFTMEAMDFTMEAMDFTMEAMSi5.445.434.494.63157164626775799399Si7Au5.505.554.244.36123 (131a)1329359325210384aThe corrected DFT value was used for potential fitting.The elastic constants of Si7Au, C11, C12, and C44 could not be predicted with equal accuracy, despite adjusting the parameter set or altering the reference structure, indicating the intrinsic limitation of the empirical model. Many studies have reported that the MEAM potential cannot well predict the elastic constants of alloys [14, 25, 30]. Therefore, C11 is usually emphasized when fitting the potential. The DFT value of 123 GPa was not directly used for parameterization owing to the nonignorable differences between the DTF-calculated and experimentally determined C11 values for Si. The MEAM potential for Si was fitted to the experimental elastic constants. Therefore, the correction (~8 GPa) to the DFT value is necessary. The error correction proposed by Ryu and Cai [25] averaged the differences between experimental and DFT results for pure Si and Au.Table5 lists the formation energies of Au introduced in Si as a substitutional or interstitial atom. The MEAM predictions agreed well with DFT results, see Ref. [14] for energy calculation. Based on Dong’s first principle study [31], a stable bridge interstitial site exists in addition to the hexagonal (Figure 1(b)) and tetrahedral (Figure 1(c)) sites. However, our DFT calculations showed that the bridge interstitial was metastable as the Au atom shifted to the tetrahedral site during relaxation. Therefore, it was not considered in the potential fitting.Table 5
Formation energy (eV) of substitutional and interstitial Au atoms in bulk Si compared to DFT values. The reference ground states for Si and Au are diamond and FCC, respectively. 64- and 216-atom systems are substitutional structures, and 65- and 217-atom systems are interstitial structures.
Point defectsDFTMEAM64/65 atoms216/217 atomsSubstitutional Au1.34a1.51a, 1.62b1.55Interstitial AuHexagonal2.28a2.37a, 2.61b2.331Tetrahedral2.32a2.39a, 2.60b2.332aCurrent work. b[31].Figure 1
Labels used to characterize the lattice distortion induced by (a) substitutional, (b) hexagonal, and (c) tetrahedral Au.
(a)(b)(c)An accurate description of the lattice distortion of the Au-doped Si crystal is required for a robust potential. The changes in the bond lengths and bond angles induced by the substitutional and interstitial Au atoms (Figure1) were examined, where n1 and n2 represent the first and second nearest neighbor atoms. The results are listed in Table 6. The distortion can be quantified by comparing the calculated values with the bond length of 2.35 Å and bond angle of 109.47° in pure Si. The lattice distortion is mainly reflected by the change in bond length l1 between Au and the first nearest neighbor Si. DFT and MEAM can predict the growth of l1, in accordance with the increased lattice constant of Si7Au.Table 6
Geometrical information corresponding to Figure1, with bond lengths in Å and angles in °. The bond length and angle are 2.35 Å and 109.47°, respectively, for pure Si.
DefectLength/angleDFTMEAMSubstitutionl12.432.47l22.352.34α109.48109.47β108.34106.96γ110.58111.87Hexagonal sitel12.452.41l22.552.53l32.322.33α62.7063.37β58.6558.31γ108.57106.89ω95.9198.67θ116.07115.23φ110.31110.92σ108.43107.36Tetrahedral sitel12.502.47l22.332.30l32.442.45α109.47109.47β68.2369.89γ114.08110.75ω111.78108.83The MEAM potential can partially describe the surface reconstruction of SiNWs without H termination. The (2×1) reconstruction on the (1 0 0) facets can be reproduced, as shown in Figures 2(a) and 2(c). However, the MEAM model cannot correctly describe the undulating surface reconstruction on the (1 1 0) facets, see the comparison between Figures 2(b) and 2(c). This unique reconstruction stems from the semiconducting nature of this material. Incorrect descriptions may lead to unpredictable deviations in the simulation results. Nevertheless, some other important surface properties can be reproduced using the MEAM potential. For example, the outermost surface shrinks toward the subsurface, as observed in another study [16].Figure 2
Unrelaxed cross-section of the ultrathin [1 0 0] SiNW (a) before and (d) after H termination. (b) Undulating reconstruction on (1 1 0) facets and (2×1) reconstruction on (1 0 0) facets predicted by DFT calculations, while (c) the MEAM potential can only predict (2×1) reconstruction. (e) and (f) Comparison of the reconstructions showing that H termination can eliminate the disparity in description on reconstruction.To eliminate the disparity in surface reconstruction, the dangling bonds were saturated with H. Figures2(e) and 2(f) and demonstrate the same surface reconstruction pattern, indicating the necessity of H termination for this study.


#### 2.2.2. Si-H Interaction

Our simulations only require H-covered SiNWs, while the interactions between H and bulk Si with interstitial H or H2 were not considered. Therefore, the MEAM parameter optimization of the Si-H binary system focused on the structural and vibrational properties of the Si-H compounds and H-covered Si surfaces. Table3 lists the MEAM parameters of the Si-H binary system. Notably, diamond was selected as the reference structure to obtain better thermal stability and vibrational characteristics.Table7 lists the structural parameters of the Si-H molecules predicted by MEAM, including the bond length, bond angle, and vibrational wave number. The Si-H bond energies in the molecules were not calculated because the reference structure of the Si-H binary system was not a dimer. Within the framework of the MEAM model, the bond energy was obtained based on the reference structure of diamond, which deviated significantly from the value obtained using dimeric SiH as the reference structure. The bond energy of dimer-Si2 predicted by MEAM was 3.87 eV, which differed significantly from the theoretical value of 3.13 eV [32]. In this circumstance, we turned the description of energy to the vibrational wave number as it measures the change of energy as a function of the bond length. Furthermore, the vibrational wave numbers predicted by MEAM agree well with the reference values. In structural properties, the MEAM estimations also agree well with the reference values, except for the bond angles θH-Si-Si in Si2H4, θH-Si-H in SiH2, and θH-Si-Si in SiH3SiH, probably owing to the inadequate description of the π-bonding.Table 7
Comparison of the structural information (bond length in Å, angle in °) and vibrational wave number (cm-1) for Si-H molecules, predicted by MEAM, experimental and theoretical methods, and other empirical calculations.
PropertyExpt./theoryTersoff [32]MEAMH2lH−H0.74a0.740.74ν4400a44004376SiH4lSi−H1.48a1.4751.49θH-Si-H109.5a109.5109.5ν2186a21002094Si2H6lSi−H1.49a, 1.50b1.481.50lSi−Si2.40a, 2.33b2.332.35θH-Si-Si110.3a, 110.1b109.4108.2ν2155a21502108Si2H4lSi−H1.48a, 1.50b1.481.52lSi−Si2.33a, 2.31b2.402.32θH-Si-H112.8a, 114.3b106.0108.7θH-Si-Si118.9a, 120.1b106.0109.3ν——2102SiHlSi−H1.51a, 1.55b1.511.56ν2042a20342145SiH2lSi−H1.51a, 1.54b1.511.52θH-Si-H92.1a, 89.8b92.3110.3ν2032a21362122SiH3lSi−H1.48a, 1.50b1.481.50θH-Si-H110.5a, 110.7b106.0110.1ν1955a20512078SiH3SiHlSi−Si2.39a, 2.38b2.372.28lSi−H (1c)1.48a, 1.51b1.481.50lSi−H (2c)1.51a, 1.54b1.511.55θH-Si-Si (1)114.2a, 114.0b109.5108.8θH-Si-Si (2)89.1a, 93.1b92.3109.2aData calculated in [32] or quoted therein. bPresent DFT calculations. c1 and 2 represent the SiH3 and SiH clusters, respectively.Regarding the SiNWs, the properties of H-covered surfaces are important. Therefore, the energy as well as the structural and vibrational properties of the H-terminated surfaces of c-Si were investigated for verification.The reconstruction pattern of the (1 0 0) surface mainly depends on the H coverage. Figure3(a)–3(c) illustrate the (2×1) monohydride, (3×1) alternate monohydride/dihydride, and (1×1) dihydriden structures produced by the MEAM potential, in line with experimental observations. As for the structural details, the Si-Si bond length on the surface of the (2×1) structure was 2.44 Å, consistent with the value of 2.4 Å [33] obtained by first-principle calculations. In addition, the bond angle of H-Si-Si 112° agrees well with other predictions 112.4° [32] and 114.7° [34]. In the (3×1) structure, the Si-Si bond length in the (2×1) part was 2.97 Å, significantly higher than the first principle prediction 2.41 Å [33]. The bond angle of H-Si-H in the (1×1) part was 112°, closed to the first principle prediction 109° [33]. In the (1×1) structure, the H-Si-H bond angle was 108°, larger than the first principle prediction 102°. Contrary to experimental observations, the asymmetric (1×1) reconstruction shown in Figure 3(c) is more stable than its symmetric counterpart, albeit close in energy [35]. In conclusion, the optimized MEAM potential provides a more accurate description of the (2×1) reconstruction. Therefore, the (1 0 0) facets of the specimens tested in the next section were modeled in the (2×1) reconstruction.Figure 3
MEAM predicted reconstructions of the (1 0 0) surface with different H coverage. (a) (2×1) monohydride. (b) Alternate (3×1) monohydride/dihydride. (c) Asymmetric (1×1) dihydride.
(a)(b)(c)The wave number of the stretch mode of (1 0 0) monohydride is 2099 cm-1, measured experimentally [36], and the (1×1) H coverage on the (1 1 1) surface gives a stretch mode wave number 2084 cm-1. The corresponding MEAM predictions of 2191 cm-1 and 2068 cm-1 agree well with the experimental values.The formation energies of H covering on the Si surfaces were considered in the parameterization, defined as follows:(5)Ef=Etot−nSiESi−nHEHnbond,where Etot is the total energy of the Si surface terminated by H, ESi and EH are the atomic average energy under the ground state, nSi and nH are the numbers of Si and H atoms in the system, and nbond is the number of Si-H bonds. The model used for calculation had a surface with 2×2 units and 3 crystalline units in thickness. Table 8 compares the formation energies calculated using MEAM and DFT.Table 8
Formation energy (eV) of H termination on Si surface.
H covered surfaceSurface energyDFTMEAM(1 0 0)-(2×1) reconstruction0.190.11(1 1 0)-monohydride0.300.27(1 1 1)-monohydride0.290.25


#### 2.2.3. H-Au Interaction

Au atoms are located within the SiNWs, and they have no complicated interactions with H atoms. Therefore, the parameterization of the H-Au binary system was conducted considering only the relationship between the potential energy of a dimeric AuH system and the atomic distance. The optimized parameters are listed in Table3, and the corresponding potential-distance curve is plotted in Figure 4. Within the range 0.45~3.7 Å, the potential energy is perfectly reproduced by MEAM.Figure 4
Potential energy as a function of atomic distance, calculated using the MEAM potential and DFT.


### 2.3. Si-H-Au Interaction

During the dynamic vibrations of SiNWs, the interactions of the H atoms with the Au-doped NWs are not complex because the H atoms are located outside the surface. The diffusion of H atoms into the NWs is beyond the scope of this study. However, when the Au atoms approach the surface, particularly in the cutoff distance from the H atoms, the force field of the ternary system works. The six parameters required to obtain the final MEAM potential of the Si-H-Au ternary system are listed in Table9. They characterize the degree of screening of the interactions between the third element and the adjacent two elements. The values were calculated using Jang’s method [30], based on the premise that the degree of screening by atom Y (in X-Y-Z) to the interaction between atoms X and Z is the average between those to the X-X and Z-Z interactions.Table 9
Parameters for the Si-H-Au ternary system.
ParameterValueCalculation formulaCmin (Si-H-Au)2.100.5CminSi‐H‐Si1/2+0.5CminAu‐H‐Au1/22Cmin (Si-Au-H)2.000.5CminSi‐Si‐H1/2+0.5CminAu‐Au‐H1/22Cmin (H-Si-Au)2.000.5CminH‐H‐Au1/2+0.5CminSi‐Si‐Au1/22Cmax (Si-H-Au)3.140.5CmaxSi‐H‐Si1/2+0.5CmaxAu‐H‐Au1/22Cmax (Si-Au-H)3.280.5CmaxSi‐Si‐H1/2+0.5CmaxAu‐Au‐H1/22Cmax (H-Si-Au)2.800.5CmaxH‐H‐Au1/2+0.5CmaxSi‐Si‐Au1/22To verify the reliability of the parameter set, the formation energies of a substitutional Au atom located at the sub-surface of the H-covered surface were calculated and compared with those obtained using DFT, see Table10. Although the MEAM parameters were not optimized targeting the DFT values, the MEAM predicted values agreed well with the reference values. Au located at the outermost surface was not considered because a previous DFT study [37] showed that the local reconstruction caused by Au was rather complex, which is beyond the description of MEAM.Table 10
Formation energies (eV) of Au located at the subsurface of the Si surfaces with H termination. The models considered here are the same as those in Table8.
SubsurfaceFormation energyDFTMEAM(1 0 0)-(2×1) reconstruction1.261.43(1 1 0)-Monohydride1.551.54(1 1 1)-Monohydride1.391.51


## 2.1. MEAM Potential for Au

Several available MEAM potentials for Au [6, 19, 22, 23] were not chosen for our work, as the parameters had to be optimized for obtaining better predictions of the fundamental properties of the Si-Au binary system. The properties referred for optimizing the Au unary MEAM potential include lattice constant, elastic constants, surface energy, and formation energy of vacancies. The MEAM values of these properties are satisfactory, as shown in Table 2. The experimental datum of the (1 0 0) surface energy was considered in the fitting. Therefore, the MEAM prediction showed that the energy of the (1 0 0) surface was much lower than that of the (1 1 0) surface. This is inconsistent with the first principle calculations but agrees with the tight-binding calculations. In addition, the MEAM potential provides higher values of the structural energy differences in face-centered cubic (FCC)/body-centered cubic (BCC) and FCC/hexagonal close-packed (HCP) structures compared to the other calculation methods, suggesting that the current MEAM potential provides a more stable FCC structure.Table 2
MEAM predictions of the lattice constanta0 (Å), elastic constants Cij (GPa), formation energy of vacancy (eV), surface energies (eV), and structural energy differences EFCC/BCC and EFCC/HCP (eV), in comparison with experimental values and higher-level calculations. The energies of vacancy formation and surface are obtained after relaxation.
PropertiesExpt.First-principleTight-bindingMEAMOthersThis worka04.07e4.068e, 4.079c—4.073e4.078C11201.63a217b, 187c—201.5b204.75C12169.67a171b, 168c—169.7b167.20C4445.44a47b, 37c—45.4b, 38d42.02B180.32a174c——179.72Ev0.89±0.04d0.55b0.49d0.91b, 0.7d0.89γ(100)1.54b1.97b1.54d1.08b, 0.95d1.55γ(110)—1.10b1.85d1.05b, 1.00d1.83γ(111)—0.92b1.48d0.90b, 0.70d1.42EFCC/BCC——0.019d0.02d0.18EFCC/HCP——0.0086d0.021d0.09a[24], 0 K values extrapolated from smooth curves, which are determined experimentally. bValues calculated in Ref. [23] or quoted therein. cDensity functional theory (DFT) calculations in this work. dValues calculated in Ref. [19] or quoted therein. eValues calculated in Ref. [25] or quoted therein.Next, we performed density functional theory (DFT) calculations to obtain the necessary reference data for developing binary and ternary MEAM potentials. The exchange correlation functional is approximated by generalized gradient approximation (GGA) revised for solids by Perdew-Burke-Ernzerhof [26], with plane wave expansion up to the cutoff energy of 300 eV determined by convergence tests. The electronic structure as a function of the positions of atomic nuclei was tackeled using the projector augmented wave method. The relaxation was repeated until the total energy reached below 10-5 eV and force below 0.01 eV/Å. The Brillouin zone was sampled with the Monkhost-Pack scheme. For calculating the cohesive energy, lattice constant, and elastic constants, 10×10×10k-points were used in the model with 1×1×1 unit and 6×6×6k-points for point defects in 3×3×3 units. The computations were completed using Vienna Ab initio Simulation Package (VASP) [27–29].Some basic physical properties of Si and Au were calculated and compared with the experimental values (given in parentheses) to verify the reliability of the above simulation settings. The lattice constant of the Si crystal was 5.439 Å (5.431 Å), the cohesive energy 4.49 eV (4.63 eV), and the elastic constants C11, C12, and C44 were 157 GPa (166 GPa), 62 GPa (64 GPa), and 75 GPa (80 GPa), respectively. The calculated and experimental values for Au are listed in Table 2.


## 2.2. MEAM Potentials of Si-Au, Si-H, and H-Au Binary Systems

### 2.2.1. Si-Au Interaction

For the Si-Au binary system, it is necessary to precisely predict the lattice distortion and changes in the elastic and energy properties of bulk Si when Au is introduced. The Si7Au compound was modeled by replacing one Si atom in the cubic unit with an Au atom. This structure is considered for parameter tuning as it represents the bulk Si doped with a high Au concentration of 12.5%. Specifically, the lattice constant, cohesive energy, and elastic constants of Si7Au and formation energies of the substitutional and interstitial Au in bulk Si were used in the parameterization. Owing to the absence of stable compound structures in the Si/Au phase diagram, b1-NaCl was chosen as the reference structure for the Si-Au binary system. Some other reference structures were considered but failed to yield predictions close to the reference values. As it is a hypothetical structure, the lattice constant and cohesive energy of b1-SiAu were adjusted slightly based on DFT data to obtain overall agreement with the reference data. Table 3 lists the MEAM parameters of the binary Si-Au system.Table 3
MEAM parameters of Si-Au, Si-H, and Au-H binary systems. The units of cohesive energyEc and nearest neighbor distance re are eV and Å, respectively. Refer to [15] for the definition of each parameter.
ParameterSi-AuSi-HH-AuReference structureb1-SiAuDiamond-SiHDimer-HAuEc5.653.492.20re1.851.5851.50α2.1332.7743.990Cmin (A-B-A)2.002.202.00Cmin (A-A-B)2.001.002.00Cmin (A-B-B)2.001.002.00Cmin (B-A-B)2.002.002.00Cmax (A-B-A)3.803.502.80Cmax (A-A-B)2.503.802.80Cmax (A-B-B)2.802.802.80Cmax (B-A-B)2.802.802.80dattrac000.18drepul04.25-0.015The performance of the potential for the Si-Au system was examined. First, the differences in the lattice constant, cohesive energy, and elastic constants between Si and Si7Au were measured. The cohesive energy was calculated as
(4)Ec=NSiESi+NAuEAu−EtotNSi+NAu,where ESi and EAu are the free energies of Si and Au atoms, respectively, NSi and NAu are the numbers of Si and Au atoms in the compound, respectively, and Etot is the total energy.The results obtained from DFT and MEAM are compared in Table4. When the concentration of Au is as high as 12.5%, the lattice distortion and energy shift can be accurately determined. The lattice constant of Si7Au increased because of the repulsion between Si and Au atoms.Table 4
Comparisons of the lattice constanta0 (Å), cohesive energy Ec (eV), and elastic constants Cij (GPa) between Si7Au and Si.
a0EcC11C12C44BDFTMEAMDFTMEAMDFTMEAMDFTMEAMDFTMEAMDFTMEAMSi5.445.434.494.63157164626775799399Si7Au5.505.554.244.36123 (131a)1329359325210384aThe corrected DFT value was used for potential fitting.The elastic constants of Si7Au, C11, C12, and C44 could not be predicted with equal accuracy, despite adjusting the parameter set or altering the reference structure, indicating the intrinsic limitation of the empirical model. Many studies have reported that the MEAM potential cannot well predict the elastic constants of alloys [14, 25, 30]. Therefore, C11 is usually emphasized when fitting the potential. The DFT value of 123 GPa was not directly used for parameterization owing to the nonignorable differences between the DTF-calculated and experimentally determined C11 values for Si. The MEAM potential for Si was fitted to the experimental elastic constants. Therefore, the correction (~8 GPa) to the DFT value is necessary. The error correction proposed by Ryu and Cai [25] averaged the differences between experimental and DFT results for pure Si and Au.Table5 lists the formation energies of Au introduced in Si as a substitutional or interstitial atom. The MEAM predictions agreed well with DFT results, see Ref. [14] for energy calculation. Based on Dong’s first principle study [31], a stable bridge interstitial site exists in addition to the hexagonal (Figure 1(b)) and tetrahedral (Figure 1(c)) sites. However, our DFT calculations showed that the bridge interstitial was metastable as the Au atom shifted to the tetrahedral site during relaxation. Therefore, it was not considered in the potential fitting.Table 5
Formation energy (eV) of substitutional and interstitial Au atoms in bulk Si compared to DFT values. The reference ground states for Si and Au are diamond and FCC, respectively. 64- and 216-atom systems are substitutional structures, and 65- and 217-atom systems are interstitial structures.
Point defectsDFTMEAM64/65 atoms216/217 atomsSubstitutional Au1.34a1.51a, 1.62b1.55Interstitial AuHexagonal2.28a2.37a, 2.61b2.331Tetrahedral2.32a2.39a, 2.60b2.332aCurrent work. b[31].Figure 1
Labels used to characterize the lattice distortion induced by (a) substitutional, (b) hexagonal, and (c) tetrahedral Au.
(a)(b)(c)An accurate description of the lattice distortion of the Au-doped Si crystal is required for a robust potential. The changes in the bond lengths and bond angles induced by the substitutional and interstitial Au atoms (Figure1) were examined, where n1 and n2 represent the first and second nearest neighbor atoms. The results are listed in Table 6. The distortion can be quantified by comparing the calculated values with the bond length of 2.35 Å and bond angle of 109.47° in pure Si. The lattice distortion is mainly reflected by the change in bond length l1 between Au and the first nearest neighbor Si. DFT and MEAM can predict the growth of l1, in accordance with the increased lattice constant of Si7Au.Table 6
Geometrical information corresponding to Figure1, with bond lengths in Å and angles in °. The bond length and angle are 2.35 Å and 109.47°, respectively, for pure Si.
DefectLength/angleDFTMEAMSubstitutionl12.432.47l22.352.34α109.48109.47β108.34106.96γ110.58111.87Hexagonal sitel12.452.41l22.552.53l32.322.33α62.7063.37β58.6558.31γ108.57106.89ω95.9198.67θ116.07115.23φ110.31110.92σ108.43107.36Tetrahedral sitel12.502.47l22.332.30l32.442.45α109.47109.47β68.2369.89γ114.08110.75ω111.78108.83The MEAM potential can partially describe the surface reconstruction of SiNWs without H termination. The (2×1) reconstruction on the (1 0 0) facets can be reproduced, as shown in Figures 2(a) and 2(c). However, the MEAM model cannot correctly describe the undulating surface reconstruction on the (1 1 0) facets, see the comparison between Figures 2(b) and 2(c). This unique reconstruction stems from the semiconducting nature of this material. Incorrect descriptions may lead to unpredictable deviations in the simulation results. Nevertheless, some other important surface properties can be reproduced using the MEAM potential. For example, the outermost surface shrinks toward the subsurface, as observed in another study [16].Figure 2
Unrelaxed cross-section of the ultrathin [1 0 0] SiNW (a) before and (d) after H termination. (b) Undulating reconstruction on (1 1 0) facets and (2×1) reconstruction on (1 0 0) facets predicted by DFT calculations, while (c) the MEAM potential can only predict (2×1) reconstruction. (e) and (f) Comparison of the reconstructions showing that H termination can eliminate the disparity in description on reconstruction.To eliminate the disparity in surface reconstruction, the dangling bonds were saturated with H. Figures2(e) and 2(f) and demonstrate the same surface reconstruction pattern, indicating the necessity of H termination for this study.


### 2.2.2. Si-H Interaction

Our simulations only require H-covered SiNWs, while the interactions between H and bulk Si with interstitial H or H2 were not considered. Therefore, the MEAM parameter optimization of the Si-H binary system focused on the structural and vibrational properties of the Si-H compounds and H-covered Si surfaces. Table3 lists the MEAM parameters of the Si-H binary system. Notably, diamond was selected as the reference structure to obtain better thermal stability and vibrational characteristics.Table7 lists the structural parameters of the Si-H molecules predicted by MEAM, including the bond length, bond angle, and vibrational wave number. The Si-H bond energies in the molecules were not calculated because the reference structure of the Si-H binary system was not a dimer. Within the framework of the MEAM model, the bond energy was obtained based on the reference structure of diamond, which deviated significantly from the value obtained using dimeric SiH as the reference structure. The bond energy of dimer-Si2 predicted by MEAM was 3.87 eV, which differed significantly from the theoretical value of 3.13 eV [32]. In this circumstance, we turned the description of energy to the vibrational wave number as it measures the change of energy as a function of the bond length. Furthermore, the vibrational wave numbers predicted by MEAM agree well with the reference values. In structural properties, the MEAM estimations also agree well with the reference values, except for the bond angles θH-Si-Si in Si2H4, θH-Si-H in SiH2, and θH-Si-Si in SiH3SiH, probably owing to the inadequate description of the π-bonding.Table 7
Comparison of the structural information (bond length in Å, angle in °) and vibrational wave number (cm-1) for Si-H molecules, predicted by MEAM, experimental and theoretical methods, and other empirical calculations.
PropertyExpt./theoryTersoff [32]MEAMH2lH−H0.74a0.740.74ν4400a44004376SiH4lSi−H1.48a1.4751.49θH-Si-H109.5a109.5109.5ν2186a21002094Si2H6lSi−H1.49a, 1.50b1.481.50lSi−Si2.40a, 2.33b2.332.35θH-Si-Si110.3a, 110.1b109.4108.2ν2155a21502108Si2H4lSi−H1.48a, 1.50b1.481.52lSi−Si2.33a, 2.31b2.402.32θH-Si-H112.8a, 114.3b106.0108.7θH-Si-Si118.9a, 120.1b106.0109.3ν——2102SiHlSi−H1.51a, 1.55b1.511.56ν2042a20342145SiH2lSi−H1.51a, 1.54b1.511.52θH-Si-H92.1a, 89.8b92.3110.3ν2032a21362122SiH3lSi−H1.48a, 1.50b1.481.50θH-Si-H110.5a, 110.7b106.0110.1ν1955a20512078SiH3SiHlSi−Si2.39a, 2.38b2.372.28lSi−H (1c)1.48a, 1.51b1.481.50lSi−H (2c)1.51a, 1.54b1.511.55θH-Si-Si (1)114.2a, 114.0b109.5108.8θH-Si-Si (2)89.1a, 93.1b92.3109.2aData calculated in [32] or quoted therein. bPresent DFT calculations. c1 and 2 represent the SiH3 and SiH clusters, respectively.Regarding the SiNWs, the properties of H-covered surfaces are important. Therefore, the energy as well as the structural and vibrational properties of the H-terminated surfaces of c-Si were investigated for verification.The reconstruction pattern of the (1 0 0) surface mainly depends on the H coverage. Figure3(a)–3(c) illustrate the (2×1) monohydride, (3×1) alternate monohydride/dihydride, and (1×1) dihydriden structures produced by the MEAM potential, in line with experimental observations. As for the structural details, the Si-Si bond length on the surface of the (2×1) structure was 2.44 Å, consistent with the value of 2.4 Å [33] obtained by first-principle calculations. In addition, the bond angle of H-Si-Si 112° agrees well with other predictions 112.4° [32] and 114.7° [34]. In the (3×1) structure, the Si-Si bond length in the (2×1) part was 2.97 Å, significantly higher than the first principle prediction 2.41 Å [33]. The bond angle of H-Si-H in the (1×1) part was 112°, closed to the first principle prediction 109° [33]. In the (1×1) structure, the H-Si-H bond angle was 108°, larger than the first principle prediction 102°. Contrary to experimental observations, the asymmetric (1×1) reconstruction shown in Figure 3(c) is more stable than its symmetric counterpart, albeit close in energy [35]. In conclusion, the optimized MEAM potential provides a more accurate description of the (2×1) reconstruction. Therefore, the (1 0 0) facets of the specimens tested in the next section were modeled in the (2×1) reconstruction.Figure 3
MEAM predicted reconstructions of the (1 0 0) surface with different H coverage. (a) (2×1) monohydride. (b) Alternate (3×1) monohydride/dihydride. (c) Asymmetric (1×1) dihydride.
(a)(b)(c)The wave number of the stretch mode of (1 0 0) monohydride is 2099 cm-1, measured experimentally [36], and the (1×1) H coverage on the (1 1 1) surface gives a stretch mode wave number 2084 cm-1. The corresponding MEAM predictions of 2191 cm-1 and 2068 cm-1 agree well with the experimental values.The formation energies of H covering on the Si surfaces were considered in the parameterization, defined as follows:(5)Ef=Etot−nSiESi−nHEHnbond,where Etot is the total energy of the Si surface terminated by H, ESi and EH are the atomic average energy under the ground state, nSi and nH are the numbers of Si and H atoms in the system, and nbond is the number of Si-H bonds. The model used for calculation had a surface with 2×2 units and 3 crystalline units in thickness. Table 8 compares the formation energies calculated using MEAM and DFT.Table 8
Formation energy (eV) of H termination on Si surface.
H covered surfaceSurface energyDFTMEAM(1 0 0)-(2×1) reconstruction0.190.11(1 1 0)-monohydride0.300.27(1 1 1)-monohydride0.290.25


### 2.2.3. H-Au Interaction

Au atoms are located within the SiNWs, and they have no complicated interactions with H atoms. Therefore, the parameterization of the H-Au binary system was conducted considering only the relationship between the potential energy of a dimeric AuH system and the atomic distance. The optimized parameters are listed in Table3, and the corresponding potential-distance curve is plotted in Figure 4. Within the range 0.45~3.7 Å, the potential energy is perfectly reproduced by MEAM.Figure 4
Potential energy as a function of atomic distance, calculated using the MEAM potential and DFT.


## 2.2.1. Si-Au Interaction

For the Si-Au binary system, it is necessary to precisely predict the lattice distortion and changes in the elastic and energy properties of bulk Si when Au is introduced. The Si7Au compound was modeled by replacing one Si atom in the cubic unit with an Au atom. This structure is considered for parameter tuning as it represents the bulk Si doped with a high Au concentration of 12.5%. Specifically, the lattice constant, cohesive energy, and elastic constants of Si7Au and formation energies of the substitutional and interstitial Au in bulk Si were used in the parameterization. Owing to the absence of stable compound structures in the Si/Au phase diagram, b1-NaCl was chosen as the reference structure for the Si-Au binary system. Some other reference structures were considered but failed to yield predictions close to the reference values. As it is a hypothetical structure, the lattice constant and cohesive energy of b1-SiAu were adjusted slightly based on DFT data to obtain overall agreement with the reference data. Table 3 lists the MEAM parameters of the binary Si-Au system.Table 3
MEAM parameters of Si-Au, Si-H, and Au-H binary systems. The units of cohesive energyEc and nearest neighbor distance re are eV and Å, respectively. Refer to [15] for the definition of each parameter.
ParameterSi-AuSi-HH-AuReference structureb1-SiAuDiamond-SiHDimer-HAuEc5.653.492.20re1.851.5851.50α2.1332.7743.990Cmin (A-B-A)2.002.202.00Cmin (A-A-B)2.001.002.00Cmin (A-B-B)2.001.002.00Cmin (B-A-B)2.002.002.00Cmax (A-B-A)3.803.502.80Cmax (A-A-B)2.503.802.80Cmax (A-B-B)2.802.802.80Cmax (B-A-B)2.802.802.80dattrac000.18drepul04.25-0.015The performance of the potential for the Si-Au system was examined. First, the differences in the lattice constant, cohesive energy, and elastic constants between Si and Si7Au were measured. The cohesive energy was calculated as
(4)Ec=NSiESi+NAuEAu−EtotNSi+NAu,where ESi and EAu are the free energies of Si and Au atoms, respectively, NSi and NAu are the numbers of Si and Au atoms in the compound, respectively, and Etot is the total energy.The results obtained from DFT and MEAM are compared in Table4. When the concentration of Au is as high as 12.5%, the lattice distortion and energy shift can be accurately determined. The lattice constant of Si7Au increased because of the repulsion between Si and Au atoms.Table 4
Comparisons of the lattice constanta0 (Å), cohesive energy Ec (eV), and elastic constants Cij (GPa) between Si7Au and Si.
a0EcC11C12C44BDFTMEAMDFTMEAMDFTMEAMDFTMEAMDFTMEAMDFTMEAMSi5.445.434.494.63157164626775799399Si7Au5.505.554.244.36123 (131a)1329359325210384aThe corrected DFT value was used for potential fitting.The elastic constants of Si7Au, C11, C12, and C44 could not be predicted with equal accuracy, despite adjusting the parameter set or altering the reference structure, indicating the intrinsic limitation of the empirical model. Many studies have reported that the MEAM potential cannot well predict the elastic constants of alloys [14, 25, 30]. Therefore, C11 is usually emphasized when fitting the potential. The DFT value of 123 GPa was not directly used for parameterization owing to the nonignorable differences between the DTF-calculated and experimentally determined C11 values for Si. The MEAM potential for Si was fitted to the experimental elastic constants. Therefore, the correction (~8 GPa) to the DFT value is necessary. The error correction proposed by Ryu and Cai [25] averaged the differences between experimental and DFT results for pure Si and Au.Table5 lists the formation energies of Au introduced in Si as a substitutional or interstitial atom. The MEAM predictions agreed well with DFT results, see Ref. [14] for energy calculation. Based on Dong’s first principle study [31], a stable bridge interstitial site exists in addition to the hexagonal (Figure 1(b)) and tetrahedral (Figure 1(c)) sites. However, our DFT calculations showed that the bridge interstitial was metastable as the Au atom shifted to the tetrahedral site during relaxation. Therefore, it was not considered in the potential fitting.Table 5
Formation energy (eV) of substitutional and interstitial Au atoms in bulk Si compared to DFT values. The reference ground states for Si and Au are diamond and FCC, respectively. 64- and 216-atom systems are substitutional structures, and 65- and 217-atom systems are interstitial structures.
Point defectsDFTMEAM64/65 atoms216/217 atomsSubstitutional Au1.34a1.51a, 1.62b1.55Interstitial AuHexagonal2.28a2.37a, 2.61b2.331Tetrahedral2.32a2.39a, 2.60b2.332aCurrent work. b[31].Figure 1
Labels used to characterize the lattice distortion induced by (a) substitutional, (b) hexagonal, and (c) tetrahedral Au.
(a)(b)(c)An accurate description of the lattice distortion of the Au-doped Si crystal is required for a robust potential. The changes in the bond lengths and bond angles induced by the substitutional and interstitial Au atoms (Figure1) were examined, where n1 and n2 represent the first and second nearest neighbor atoms. The results are listed in Table 6. The distortion can be quantified by comparing the calculated values with the bond length of 2.35 Å and bond angle of 109.47° in pure Si. The lattice distortion is mainly reflected by the change in bond length l1 between Au and the first nearest neighbor Si. DFT and MEAM can predict the growth of l1, in accordance with the increased lattice constant of Si7Au.Table 6
Geometrical information corresponding to Figure1, with bond lengths in Å and angles in °. The bond length and angle are 2.35 Å and 109.47°, respectively, for pure Si.
DefectLength/angleDFTMEAMSubstitutionl12.432.47l22.352.34α109.48109.47β108.34106.96γ110.58111.87Hexagonal sitel12.452.41l22.552.53l32.322.33α62.7063.37β58.6558.31γ108.57106.89ω95.9198.67θ116.07115.23φ110.31110.92σ108.43107.36Tetrahedral sitel12.502.47l22.332.30l32.442.45α109.47109.47β68.2369.89γ114.08110.75ω111.78108.83The MEAM potential can partially describe the surface reconstruction of SiNWs without H termination. The (2×1) reconstruction on the (1 0 0) facets can be reproduced, as shown in Figures 2(a) and 2(c). However, the MEAM model cannot correctly describe the undulating surface reconstruction on the (1 1 0) facets, see the comparison between Figures 2(b) and 2(c). This unique reconstruction stems from the semiconducting nature of this material. Incorrect descriptions may lead to unpredictable deviations in the simulation results. Nevertheless, some other important surface properties can be reproduced using the MEAM potential. For example, the outermost surface shrinks toward the subsurface, as observed in another study [16].Figure 2
Unrelaxed cross-section of the ultrathin [1 0 0] SiNW (a) before and (d) after H termination. (b) Undulating reconstruction on (1 1 0) facets and (2×1) reconstruction on (1 0 0) facets predicted by DFT calculations, while (c) the MEAM potential can only predict (2×1) reconstruction. (e) and (f) Comparison of the reconstructions showing that H termination can eliminate the disparity in description on reconstruction.To eliminate the disparity in surface reconstruction, the dangling bonds were saturated with H. Figures2(e) and 2(f) and demonstrate the same surface reconstruction pattern, indicating the necessity of H termination for this study.


## 2.2.2. Si-H Interaction

Our simulations only require H-covered SiNWs, while the interactions between H and bulk Si with interstitial H or H2 were not considered. Therefore, the MEAM parameter optimization of the Si-H binary system focused on the structural and vibrational properties of the Si-H compounds and H-covered Si surfaces. Table3 lists the MEAM parameters of the Si-H binary system. Notably, diamond was selected as the reference structure to obtain better thermal stability and vibrational characteristics.Table7 lists the structural parameters of the Si-H molecules predicted by MEAM, including the bond length, bond angle, and vibrational wave number. The Si-H bond energies in the molecules were not calculated because the reference structure of the Si-H binary system was not a dimer. Within the framework of the MEAM model, the bond energy was obtained based on the reference structure of diamond, which deviated significantly from the value obtained using dimeric SiH as the reference structure. The bond energy of dimer-Si2 predicted by MEAM was 3.87 eV, which differed significantly from the theoretical value of 3.13 eV [32]. In this circumstance, we turned the description of energy to the vibrational wave number as it measures the change of energy as a function of the bond length. Furthermore, the vibrational wave numbers predicted by MEAM agree well with the reference values. In structural properties, the MEAM estimations also agree well with the reference values, except for the bond angles θH-Si-Si in Si2H4, θH-Si-H in SiH2, and θH-Si-Si in SiH3SiH, probably owing to the inadequate description of the π-bonding.Table 7
Comparison of the structural information (bond length in Å, angle in °) and vibrational wave number (cm-1) for Si-H molecules, predicted by MEAM, experimental and theoretical methods, and other empirical calculations.
PropertyExpt./theoryTersoff [32]MEAMH2lH−H0.74a0.740.74ν4400a44004376SiH4lSi−H1.48a1.4751.49θH-Si-H109.5a109.5109.5ν2186a21002094Si2H6lSi−H1.49a, 1.50b1.481.50lSi−Si2.40a, 2.33b2.332.35θH-Si-Si110.3a, 110.1b109.4108.2ν2155a21502108Si2H4lSi−H1.48a, 1.50b1.481.52lSi−Si2.33a, 2.31b2.402.32θH-Si-H112.8a, 114.3b106.0108.7θH-Si-Si118.9a, 120.1b106.0109.3ν——2102SiHlSi−H1.51a, 1.55b1.511.56ν2042a20342145SiH2lSi−H1.51a, 1.54b1.511.52θH-Si-H92.1a, 89.8b92.3110.3ν2032a21362122SiH3lSi−H1.48a, 1.50b1.481.50θH-Si-H110.5a, 110.7b106.0110.1ν1955a20512078SiH3SiHlSi−Si2.39a, 2.38b2.372.28lSi−H (1c)1.48a, 1.51b1.481.50lSi−H (2c)1.51a, 1.54b1.511.55θH-Si-Si (1)114.2a, 114.0b109.5108.8θH-Si-Si (2)89.1a, 93.1b92.3109.2aData calculated in [32] or quoted therein. bPresent DFT calculations. c1 and 2 represent the SiH3 and SiH clusters, respectively.Regarding the SiNWs, the properties of H-covered surfaces are important. Therefore, the energy as well as the structural and vibrational properties of the H-terminated surfaces of c-Si were investigated for verification.The reconstruction pattern of the (1 0 0) surface mainly depends on the H coverage. Figure3(a)–3(c) illustrate the (2×1) monohydride, (3×1) alternate monohydride/dihydride, and (1×1) dihydriden structures produced by the MEAM potential, in line with experimental observations. As for the structural details, the Si-Si bond length on the surface of the (2×1) structure was 2.44 Å, consistent with the value of 2.4 Å [33] obtained by first-principle calculations. In addition, the bond angle of H-Si-Si 112° agrees well with other predictions 112.4° [32] and 114.7° [34]. In the (3×1) structure, the Si-Si bond length in the (2×1) part was 2.97 Å, significantly higher than the first principle prediction 2.41 Å [33]. The bond angle of H-Si-H in the (1×1) part was 112°, closed to the first principle prediction 109° [33]. In the (1×1) structure, the H-Si-H bond angle was 108°, larger than the first principle prediction 102°. Contrary to experimental observations, the asymmetric (1×1) reconstruction shown in Figure 3(c) is more stable than its symmetric counterpart, albeit close in energy [35]. In conclusion, the optimized MEAM potential provides a more accurate description of the (2×1) reconstruction. Therefore, the (1 0 0) facets of the specimens tested in the next section were modeled in the (2×1) reconstruction.Figure 3
MEAM predicted reconstructions of the (1 0 0) surface with different H coverage. (a) (2×1) monohydride. (b) Alternate (3×1) monohydride/dihydride. (c) Asymmetric (1×1) dihydride.
(a)(b)(c)The wave number of the stretch mode of (1 0 0) monohydride is 2099 cm-1, measured experimentally [36], and the (1×1) H coverage on the (1 1 1) surface gives a stretch mode wave number 2084 cm-1. The corresponding MEAM predictions of 2191 cm-1 and 2068 cm-1 agree well with the experimental values.The formation energies of H covering on the Si surfaces were considered in the parameterization, defined as follows:(5)Ef=Etot−nSiESi−nHEHnbond,where Etot is the total energy of the Si surface terminated by H, ESi and EH are the atomic average energy under the ground state, nSi and nH are the numbers of Si and H atoms in the system, and nbond is the number of Si-H bonds. The model used for calculation had a surface with 2×2 units and 3 crystalline units in thickness. Table 8 compares the formation energies calculated using MEAM and DFT.Table 8
Formation energy (eV) of H termination on Si surface.
H covered surfaceSurface energyDFTMEAM(1 0 0)-(2×1) reconstruction0.190.11(1 1 0)-monohydride0.300.27(1 1 1)-monohydride0.290.25


## 2.2.3. H-Au Interaction

Au atoms are located within the SiNWs, and they have no complicated interactions with H atoms. Therefore, the parameterization of the H-Au binary system was conducted considering only the relationship between the potential energy of a dimeric AuH system and the atomic distance. The optimized parameters are listed in Table3, and the corresponding potential-distance curve is plotted in Figure 4. Within the range 0.45~3.7 Å, the potential energy is perfectly reproduced by MEAM.Figure 4
Potential energy as a function of atomic distance, calculated using the MEAM potential and DFT.


## 2.3. Si-H-Au Interaction

During the dynamic vibrations of SiNWs, the interactions of the H atoms with the Au-doped NWs are not complex because the H atoms are located outside the surface. The diffusion of H atoms into the NWs is beyond the scope of this study. However, when the Au atoms approach the surface, particularly in the cutoff distance from the H atoms, the force field of the ternary system works. The six parameters required to obtain the final MEAM potential of the Si-H-Au ternary system are listed in Table9. They characterize the degree of screening of the interactions between the third element and the adjacent two elements. The values were calculated using Jang’s method [30], based on the premise that the degree of screening by atom Y (in X-Y-Z) to the interaction between atoms X and Z is the average between those to the X-X and Z-Z interactions.Table 9
Parameters for the Si-H-Au ternary system.
ParameterValueCalculation formulaCmin (Si-H-Au)2.100.5CminSi‐H‐Si1/2+0.5CminAu‐H‐Au1/22Cmin (Si-Au-H)2.000.5CminSi‐Si‐H1/2+0.5CminAu‐Au‐H1/22Cmin (H-Si-Au)2.000.5CminH‐H‐Au1/2+0.5CminSi‐Si‐Au1/22Cmax (Si-H-Au)3.140.5CmaxSi‐H‐Si1/2+0.5CmaxAu‐H‐Au1/22Cmax (Si-Au-H)3.280.5CmaxSi‐Si‐H1/2+0.5CmaxAu‐Au‐H1/22Cmax (H-Si-Au)2.800.5CmaxH‐H‐Au1/2+0.5CmaxSi‐Si‐Au1/22To verify the reliability of the parameter set, the formation energies of a substitutional Au atom located at the sub-surface of the H-covered surface were calculated and compared with those obtained using DFT, see Table10. Although the MEAM parameters were not optimized targeting the DFT values, the MEAM predicted values agreed well with the reference values. Au located at the outermost surface was not considered because a previous DFT study [37] showed that the local reconstruction caused by Au was rather complex, which is beyond the description of MEAM.Table 10
Formation energies (eV) of Au located at the subsurface of the Si surfaces with H termination. The models considered here are the same as those in Table8.
SubsurfaceFormation energyDFTMEAM(1 0 0)-(2×1) reconstruction1.261.43(1 1 0)-Monohydride1.551.54(1 1 1)-Monohydride1.391.51


## 3. Vibrational Properties of H-Covered SiNWs Doped with Au

### 3.1. Modeling and Vibrational Test

After consulting relevant literature, we can only determine the cross-sectional shape of [1 0 0] and [1 1 0] SiNWs. The cross-section of the ultrathin [1 0 0] SiNW is octagonal with four (1 0 0) facets and four (1 1 0) facets [38], as shown in Figure 5(b). The (1 0 0) facets are in the (2×1) reconstruction. The cross-section of the ultrathin [1 1 0] SiNW is hexagonal with two (1 0 0) facets and four (1 1 1) facets [38], see Figure 5(c). Au atoms randomly replace Si atoms according to the doping concentration. Other forms of doping were not considered owing to the absence of the doping profile of Si-Au compounds. According to Liang’s first principle study, Au is likely to be located at the substitutional sites [39]. As Au and Si are not congeners in the same group, when Au is located on the surface of the NWs, it would generate complex surface reconstructions that cannot be described using the MEAM potential. Therefore, the doping region is enclosed within the subsurface. Actually first-principles calculations [39, 40] and experiments [41] suggested that the core region of SiNWs is energetically favorable for Au atoms.Figure 5
(a) Schematic of a clamped-clamped SiNW used for vibrational tests. A sinusoidal velocity field introduces the initial excitation. (b) Cross-section of the [1 0 0] SiNW. (c) Cross-section of the [1 1 0] SiNW. Si, Au, and H atoms are colored blue, orange, and yellow, respectively.
(a)(b)(c)The model size was determined according to the available computing resources. Because of the complexity of the MEAM scheme and the large cutoff distance (6 Å), it is considerably more expensive than SW, Tersoff, and EAM models. Besides, the investigation on the influence of Au on the quality factor of SiNWs requires at least 5 ns. To address the above concerns, the [1 0 0]-SiNW for vibrational tests was 3.3 nm in diameter and 18.5 nm in length, containing 9282 atoms (1496 H atoms); the [1 1 0]-SiNW was 3.4 nm in diameter and 13.8 nm in length, containing 7200 atoms (1152 H atoms). To compare the vibrational characteristics of the H-covered and bare NWs, the bare [1 0 0]-SiNWs were obtained by removing the H atoms of the H-covered specimens. Another reason for not preparing larger specimens is that the size effects on the vibrational characteristics have already been studied [12].The SiNW was first relaxed to a minimum energy state using the conjugate gradient method. Then, the SiNW was equilibrated at 10 K using Nose-Hoover thermostat (NVT ensemble) for 100 ps so that the NW could further release the residual stress and reach thermodynamic equilibrium. After static and dynamic relaxations, one lattice unit at both ends of the SiNW was fixed, and a sinusoidal velocity fieldvx=asinkx0≤x≤L was applied to the movable part, as shown in Figure 5(a). The amplitude a was 0.8 Å/ps, small enough to prevent nonlinear vibration. In the last step, the SiNW vibrated under the constant energy ensemble (NVE) for 5 ~ 10 ns. The length of time t depended on the period of beat phenomenon, Tb. If Tb was higher than t, then the vibration continued and exceeded one beat period. The coordinates of the center of mass (COM) and energy during vibration were recorded. The computations were accomplished through the large-scale atomic/molecular massively parallel simulator (LAMMPS) [42], using 2 nodes with 16 cores. Each test took 7 ~ 15 hours, depending on the time span of the vibration.


### 3.2. Natural Frequency Decrease Caused by the Impurity Au

Fast Fourier transform (FFT) was applied to the time history of the potential energy or root mean square (RMS) of the COM to obtain the first-order natural frequency. The RMS of the COM was considered for later analysis ofQ. We considered a [1 0 0] sample in the 3% group as an example to illustrate signal processing. Figure 6(a) plots the time history of the potential energy and the magnified view of the periodic vibration. The period is constant although the oscillation amplitude keeps changing. FFT provides one-sided power spectra, as shown in Figure 6(c). Although the FFT magnitude is different between these two signals, the spectrum distribution is consistent, as the two peaks in each curve correspond to the same frequencies. This indicates that the RMS of the COM can replace the energy for frequency analysis.Figure 6
Processing of the vibrational information. A [1 0 0] sample in the 3% group was taken as an example. It was initially excited by the stimuli alongy and z as a superposition. Time histories of (a) potential energy, (b) yCOM and zCOM, and (d) RMS of COM exhibit beat phenomena, reflected by the two adjacent peaks in the spectrum distribution obtained by FFT (c). The inset gives the FFT spectrum of the new elementary directions y′ and z′, along which the beat phenomena disappear. (d) Attenuation of the upper envelope corresponding to the dissipation of mechanical energy.The effect of Au on the resonant frequency was investigated considering two aspects, concentration and distribution. To this end, five samples with different impurity distribution were prepared for each concentration (≤7%). The natural frequency as a function of concentration is plotted in Figure7, which exhibits an evident linear downward trend from 0 to 7%. The error bar corresponding to each concentration is negligible compared to the frequency value, indicating that the impurity distribution has a negligible influence on the natural frequency.Figure 7
Resonant frequency of the double-clamped SiNWs as a function of doping concentration. The error bar corresponding to each concentration was sampled from five specimens with different impurity distributions.The two adjacent peaks in the spectrum indicate the vibrations in two orthogonal directions [43], as plotted in Figure 6(c), which are denoted as the elementary directions. Synthetizing the two vibrations results in the beat phenomenon [44]. Symmetry breaking in NWs leads to the deflection of the elementary directions [13]. Obviously, Au doping caused lattice asymmetry owing to the significant difference in the atomic masses of Au and Si, theoretically leading to the deflection of the elementary directions. However, the simulation results for [1 1 0] SiNWs showed that the elementary directions remained unchanged (i.e., the y and z orientations in Figure 5(c)) or deflect at a negligible angle, irrespective of different doping concentrations, distributions, or initial excitation directions. These aspects will be further explained in the next part, with this part focusing on the frequency. The frequencies along the elementary directions depend on the symmetry. The frequencies along the elementary directions y and z in the [1 0 0] SiNWs are close (see Figure 6(c)) even though the elementary axes deflect significantly, whereas those in the [1 1 0] SiNWs differ considerably. Therefore, only the values of the principal peaks in the FFT curve of the [1 0 0] SiNWs are plotted in Figure 7.The effects of H termination on the frequency were investigated. The trend for bare [1 0 0] samples is linear, similar to the H-covered counterparts, as shown in Figure7. H termination causes a slight decrement in frequency but does not affect the gradient. These results suggest the advantages of MD simulations involving the NWs with complex cross-sections over theoretical calculations. Based on the Euler-Bernoulli beam theory (EBBT), He and Lilley [10] developed a new formula to calculate the natural frequency of a double-clamped beam considering surface stress, i.e., f0=fn0.025η+1 and η=2τDL2/EI∗. fn, τ, D, L, and (EI)∗ are the frequency obtained using EBBT, surface stress, diameter, length, and effective bending moment, respectively. The equation is unsuitable for the H-covered SiNWs, because of not only the simplification of the formula itself but also the approximation of the diameter, effective bending inertia, and surface stress that makes it difficult to capture the small differences accurately.


### 3.3. Beat Phenomenon

The beat phenomenon stems from the superposition of two harmonic vibrations with close frequencies. For the convenience of theoretical analysis, the amplitude attenuation was ignored. For two harmonic vibrations described byUi=AisinωitDi=Aisinwit,, where ωi=2πfi, i=y, z, the total external energy of the NW with mass M can be expressed as [13]
(6)Ee=∑i=1,212MAi2ωi2sin2ωit=Aetot‐Aetotcosω1−ω2tcosω1+ω2t+MA12ω12−MA22ω22sinω1−ω2tsinω1+ω2t≈Aetot‐Aetotcosω1−ω2tcosω1+ω2t,where Aetot=MA12ω12+MA22ω22/4 is the amplitude of the total external energy. The approximation in Equation (6) is based on the premise ωy≈ωz and sin[(ωy-ωz)t] approaching 0. Equation Equation (6) presents a vibration with the oscillating amplitude at a low frequency, following Aetotcosωy−ωzt.The vibration decomposes into two elementary directions when the initial stimulus is along an arbitrary direction in the cross-section, resulting in a beat phenomenon [13]. Different impurity distributions create different symmetries determined by the mass distribution rather than the geometry. In this way, the elementary coordinate system changes with the distribution. Therefore, different distributions can lead to different beat phenomena because the elementary directions change. However, the theoretical possibility explained above is inconsistent with the simulation results for [1 1 0] SiNWs, the elementary directions of which remained unchanged or deflected at negligible angles. It is interesting to examine the robustness of the elementary directions. According to the model pattern used in this work, several possible factors could cause the deflection of the elementary directions, including H termination, and the impurity concentration and distribution.First, the effects of H termination were investigated. One (Figure8(a)) of the five [1 1 0] samples in the 7% group was tested after removing the H atoms. The stimulus direction was anticlockwise at π/4 from y. The time history of the COM along z and the corresponding FFT is shown in Figure 8(b). The beat phenomenon can be observed, which is characterized by the auxiliary peak, suggesting the deflection of the elementary directions. The new elementary coordinate system y′-z′ could be obtained using the following transformation matrix:
(7)y′z′=cosθ−sinθsinθcosθyz,where θ is the deflection angle between y′-z′ and y-z, and it was obtained by trial and error until no beat phenomenon appeared along y′ and z′. θ was 0.0095π in this case, obviously a negligible angle. The vibration in the transformed coordinate system is shown in Figure 8(c), without the beat phenomenon. It is worth mentioning that the H-covered counterpart of this model showed no beat phenomenon, indicating that H coverage can improve the robustness of the elementary directions to a certain extent.Figure 8
Effects of H termination (a)~(c), radial (d)~(f), and axial impurity segregation (g, f) on the deflection of the elementary directions. (b) Appearance of the beat phenomenon after removing the H atoms on the surface of a sample in the 7% group (a) and disappearance after rotating the coordinate system by 0.0095π (c). Deteriorated radial symmetry (d) resulting in the beat phenomenon (e) and its disappearance after rotating the coordinate system by 0.0064π (f). Axial symmetry breaking (g) causing no deflection of the elementary directions (h).Although the elementary directions proved robust against the impurity concentration and distribution, different beat periods for the five samples in each group with different concentrations were obtained. This is because once the elementary direction changes, the frequency along this direction changes accordingly. The beat angular frequency is (ωy-ωz). A slight change in ωy or ωz would result in a different beat frequency.The beat frequency can be calculated as follows considering an NW with a fixed cross-section and length [43].
(8)ωy−ωz∝Iy‐Iz,where I is the moment of inertia. However, this equation could not be used to estimate the beat frequency due to the randomness in the distribution of Au atoms. Therefore, symmetry breaking is difficult to quantify. Hence, this study applied an alternative method to characterize the degree of symmetry breaking, i.e., the enveloping cylindrical surface (sharing the same axis with the NW) of the axial mass distribution curve. The radius of the enveloping cylindrical surface in the bare NW (Figure 8(a)) is 1.95 Å and 1.93 Å for the H-covered counterpart. The slightly lower value indicates better symmetry.Then, the effects of the symmetry of the mass distribution were investigated, considering more extreme situations, i.e., when the impurity atoms were radially and axially segregated. Figure8(d) presents a SiNW with Au segregated in half of the region, which was constructed based on the H-covered counterpart in Figure 8(a) obtained by replacing the Au atoms in the other half region with Si atoms. The radius of the enveloping cylindrical surface increased to 2.34 Å, indicating deteriorated symmetry. The excitation direction was the same as above. As observed in Figure 8(e), the model shows a slight beat phenomenon. The auxiliary peak rotates y-z by 0.0064π to obtain the new elementary coordinate system y′-z′. The transformed vibrational signal and spectrum are plotted in Figure 8(f). The rotation angle is also negligible, indicating that radial symmetry breaking slightly affects the deflection of the elementary directions. In comparison, axial symmetry breaking will not cause the deflection of the elementary directions. The model shown in Figure 8(g) is the H-covered counterpart of Figure 9(a) with only one-third of the doped region in the axial direction. The result (Figure 8(h)) shows no beat phenomenon, suggesting that the elementary directions remain unchanged.Figure 9
(a) Quality factor of [1 0 0] and [1 1 0] and (b) beat periodTb SiNWs as a function of concentration.
(a)(b)The above investigations lead to the conclusion that the robustness of the elementary directions in the [1 1 0] SiNWs is hardly affected by impurity Au atoms. However, from a rigorous point of view, the robustness may also be related to the cross-section pattern and boundary conditions (e.g., clamped on one side) or other factors associated with impurity Au atoms, which is a topic for future work. Furthermore, the reason that [1 0 0] SiNWs do not possess this characteristic is also an interesting topic. Actually, for most Au-doped [1 0 0] samples, we failed to find the new elementary axesy′-z′ using Equation (7), irrespective of how we adjust θ, and the FFT curves for y′ and z′ always exhibited an auxiliary peak. Thus, we can preliminarily infer that this behavior is closely related to the symmetry difference and the original elementary directions.This section explains that the unintentional impurity Au atoms will not cause serious deflection of the elementary axes in the [1 1 0] SiNWs. When applied in ultrasensitive resonators, the symmetry breaking induced by the impurity distribution will not produce a large frequency shift. It will not cause disturbance during signal testing. According to the advanced mass detection technique developed by Gil Santos et al. [45], the position and mass of the adsorbate can be obtained by the relative shift of frequencies in two orthogonal directions. Effective restraints on extra frequency shifts caused by the deflection of elementary directions are required. [1 1 0] SiNWs can meet this need because their elementary axes are robust.


### 3.4. Quality FactorQ

Q quantifies the energy dissipation rate of a resonator, defined as the ratio between the total external energy and the average energy loss in one cycle at the resonant frequency. MD can output the total energy, including the external energy (or mechanical energy) and internal energy (or thermal energy). It is inconvenient to extract the external energy to obtain Q; thus, we used an alternative method, i.e., the RMS of the COM [46]. The amplitude of the RMS decreases following the exponential law Ae-ξωt, where ξ is the damping ratio, and ω =2πf0 is the angular frequency, with f0 obtained from FFT. Moreover, Q is equal to 1/2ξ.TheQ values for [1 0 0] and [1 1 0] SiNWs are plotted in Figure 9(a) as a function of concentration. The graph shows an overall downward trend, accompanied by the enlarged error bar. Besides, the impurity distribution considerably magnifies the scattering of Q. The scattering of Q is amplified as the concentration increases. In this case, we cannot safely capture the varying trend of Q as a function of the concentration due to the small sample size. In addition to external damping, atomic friction is also an important factor for energy dissipation [47]. According to Islam et al. [4], impurity atoms can promote friction, which results in decreased Q. Therefore, the downward trend is logically correct.The strong dependence of the beat period,Tb=2π/ωy−ωz, of the [1 1 0] SiNWs on the concentration is shown in Figure 9(b). The longer the beat period, the larger is the deflection of the elementary coordinate system. Although the elementary coordinate system is robust, a slight deflection can lead to an obvious beat phenomenon. From Figures 9(a) and 9(b), one can infer that the deflection of the elementary axes accelerates the vibrational energy dissipation. This feature was not observed in the [1 0 0] SiNWs, probably because of the same geometric symmetry along y and z.More computational resources are required to reveal the relationship betweenQ and the doping concentration than for the other properties, only to obtain the downward trend and increased scattering. Nevertheless, these findings have engineering significance. For example, high Au concentrations should be avoided for standardizing SiNW-based resonators.


## 3.1. Modeling and Vibrational Test

After consulting relevant literature, we can only determine the cross-sectional shape of [1 0 0] and [1 1 0] SiNWs. The cross-section of the ultrathin [1 0 0] SiNW is octagonal with four (1 0 0) facets and four (1 1 0) facets [38], as shown in Figure 5(b). The (1 0 0) facets are in the (2×1) reconstruction. The cross-section of the ultrathin [1 1 0] SiNW is hexagonal with two (1 0 0) facets and four (1 1 1) facets [38], see Figure 5(c). Au atoms randomly replace Si atoms according to the doping concentration. Other forms of doping were not considered owing to the absence of the doping profile of Si-Au compounds. According to Liang’s first principle study, Au is likely to be located at the substitutional sites [39]. As Au and Si are not congeners in the same group, when Au is located on the surface of the NWs, it would generate complex surface reconstructions that cannot be described using the MEAM potential. Therefore, the doping region is enclosed within the subsurface. Actually first-principles calculations [39, 40] and experiments [41] suggested that the core region of SiNWs is energetically favorable for Au atoms.Figure 5
(a) Schematic of a clamped-clamped SiNW used for vibrational tests. A sinusoidal velocity field introduces the initial excitation. (b) Cross-section of the [1 0 0] SiNW. (c) Cross-section of the [1 1 0] SiNW. Si, Au, and H atoms are colored blue, orange, and yellow, respectively.
(a)(b)(c)The model size was determined according to the available computing resources. Because of the complexity of the MEAM scheme and the large cutoff distance (6 Å), it is considerably more expensive than SW, Tersoff, and EAM models. Besides, the investigation on the influence of Au on the quality factor of SiNWs requires at least 5 ns. To address the above concerns, the [1 0 0]-SiNW for vibrational tests was 3.3 nm in diameter and 18.5 nm in length, containing 9282 atoms (1496 H atoms); the [1 1 0]-SiNW was 3.4 nm in diameter and 13.8 nm in length, containing 7200 atoms (1152 H atoms). To compare the vibrational characteristics of the H-covered and bare NWs, the bare [1 0 0]-SiNWs were obtained by removing the H atoms of the H-covered specimens. Another reason for not preparing larger specimens is that the size effects on the vibrational characteristics have already been studied [12].The SiNW was first relaxed to a minimum energy state using the conjugate gradient method. Then, the SiNW was equilibrated at 10 K using Nose-Hoover thermostat (NVT ensemble) for 100 ps so that the NW could further release the residual stress and reach thermodynamic equilibrium. After static and dynamic relaxations, one lattice unit at both ends of the SiNW was fixed, and a sinusoidal velocity fieldvx=asinkx0≤x≤L was applied to the movable part, as shown in Figure 5(a). The amplitude a was 0.8 Å/ps, small enough to prevent nonlinear vibration. In the last step, the SiNW vibrated under the constant energy ensemble (NVE) for 5 ~ 10 ns. The length of time t depended on the period of beat phenomenon, Tb. If Tb was higher than t, then the vibration continued and exceeded one beat period. The coordinates of the center of mass (COM) and energy during vibration were recorded. The computations were accomplished through the large-scale atomic/molecular massively parallel simulator (LAMMPS) [42], using 2 nodes with 16 cores. Each test took 7 ~ 15 hours, depending on the time span of the vibration.


## 3.2. Natural Frequency Decrease Caused by the Impurity Au

Fast Fourier transform (FFT) was applied to the time history of the potential energy or root mean square (RMS) of the COM to obtain the first-order natural frequency. The RMS of the COM was considered for later analysis ofQ. We considered a [1 0 0] sample in the 3% group as an example to illustrate signal processing. Figure 6(a) plots the time history of the potential energy and the magnified view of the periodic vibration. The period is constant although the oscillation amplitude keeps changing. FFT provides one-sided power spectra, as shown in Figure 6(c). Although the FFT magnitude is different between these two signals, the spectrum distribution is consistent, as the two peaks in each curve correspond to the same frequencies. This indicates that the RMS of the COM can replace the energy for frequency analysis.Figure 6
Processing of the vibrational information. A [1 0 0] sample in the 3% group was taken as an example. It was initially excited by the stimuli alongy and z as a superposition. Time histories of (a) potential energy, (b) yCOM and zCOM, and (d) RMS of COM exhibit beat phenomena, reflected by the two adjacent peaks in the spectrum distribution obtained by FFT (c). The inset gives the FFT spectrum of the new elementary directions y′ and z′, along which the beat phenomena disappear. (d) Attenuation of the upper envelope corresponding to the dissipation of mechanical energy.The effect of Au on the resonant frequency was investigated considering two aspects, concentration and distribution. To this end, five samples with different impurity distribution were prepared for each concentration (≤7%). The natural frequency as a function of concentration is plotted in Figure7, which exhibits an evident linear downward trend from 0 to 7%. The error bar corresponding to each concentration is negligible compared to the frequency value, indicating that the impurity distribution has a negligible influence on the natural frequency.Figure 7
Resonant frequency of the double-clamped SiNWs as a function of doping concentration. The error bar corresponding to each concentration was sampled from five specimens with different impurity distributions.The two adjacent peaks in the spectrum indicate the vibrations in two orthogonal directions [43], as plotted in Figure 6(c), which are denoted as the elementary directions. Synthetizing the two vibrations results in the beat phenomenon [44]. Symmetry breaking in NWs leads to the deflection of the elementary directions [13]. Obviously, Au doping caused lattice asymmetry owing to the significant difference in the atomic masses of Au and Si, theoretically leading to the deflection of the elementary directions. However, the simulation results for [1 1 0] SiNWs showed that the elementary directions remained unchanged (i.e., the y and z orientations in Figure 5(c)) or deflect at a negligible angle, irrespective of different doping concentrations, distributions, or initial excitation directions. These aspects will be further explained in the next part, with this part focusing on the frequency. The frequencies along the elementary directions depend on the symmetry. The frequencies along the elementary directions y and z in the [1 0 0] SiNWs are close (see Figure 6(c)) even though the elementary axes deflect significantly, whereas those in the [1 1 0] SiNWs differ considerably. Therefore, only the values of the principal peaks in the FFT curve of the [1 0 0] SiNWs are plotted in Figure 7.The effects of H termination on the frequency were investigated. The trend for bare [1 0 0] samples is linear, similar to the H-covered counterparts, as shown in Figure7. H termination causes a slight decrement in frequency but does not affect the gradient. These results suggest the advantages of MD simulations involving the NWs with complex cross-sections over theoretical calculations. Based on the Euler-Bernoulli beam theory (EBBT), He and Lilley [10] developed a new formula to calculate the natural frequency of a double-clamped beam considering surface stress, i.e., f0=fn0.025η+1 and η=2τDL2/EI∗. fn, τ, D, L, and (EI)∗ are the frequency obtained using EBBT, surface stress, diameter, length, and effective bending moment, respectively. The equation is unsuitable for the H-covered SiNWs, because of not only the simplification of the formula itself but also the approximation of the diameter, effective bending inertia, and surface stress that makes it difficult to capture the small differences accurately.


## 3.3. Beat Phenomenon

The beat phenomenon stems from the superposition of two harmonic vibrations with close frequencies. For the convenience of theoretical analysis, the amplitude attenuation was ignored. For two harmonic vibrations described byUi=AisinωitDi=Aisinwit,, where ωi=2πfi, i=y, z, the total external energy of the NW with mass M can be expressed as [13]
(6)Ee=∑i=1,212MAi2ωi2sin2ωit=Aetot‐Aetotcosω1−ω2tcosω1+ω2t+MA12ω12−MA22ω22sinω1−ω2tsinω1+ω2t≈Aetot‐Aetotcosω1−ω2tcosω1+ω2t,where Aetot=MA12ω12+MA22ω22/4 is the amplitude of the total external energy. The approximation in Equation (6) is based on the premise ωy≈ωz and sin[(ωy-ωz)t] approaching 0. Equation Equation (6) presents a vibration with the oscillating amplitude at a low frequency, following Aetotcosωy−ωzt.The vibration decomposes into two elementary directions when the initial stimulus is along an arbitrary direction in the cross-section, resulting in a beat phenomenon [13]. Different impurity distributions create different symmetries determined by the mass distribution rather than the geometry. In this way, the elementary coordinate system changes with the distribution. Therefore, different distributions can lead to different beat phenomena because the elementary directions change. However, the theoretical possibility explained above is inconsistent with the simulation results for [1 1 0] SiNWs, the elementary directions of which remained unchanged or deflected at negligible angles. It is interesting to examine the robustness of the elementary directions. According to the model pattern used in this work, several possible factors could cause the deflection of the elementary directions, including H termination, and the impurity concentration and distribution.First, the effects of H termination were investigated. One (Figure8(a)) of the five [1 1 0] samples in the 7% group was tested after removing the H atoms. The stimulus direction was anticlockwise at π/4 from y. The time history of the COM along z and the corresponding FFT is shown in Figure 8(b). The beat phenomenon can be observed, which is characterized by the auxiliary peak, suggesting the deflection of the elementary directions. The new elementary coordinate system y′-z′ could be obtained using the following transformation matrix:
(7)y′z′=cosθ−sinθsinθcosθyz,where θ is the deflection angle between y′-z′ and y-z, and it was obtained by trial and error until no beat phenomenon appeared along y′ and z′. θ was 0.0095π in this case, obviously a negligible angle. The vibration in the transformed coordinate system is shown in Figure 8(c), without the beat phenomenon. It is worth mentioning that the H-covered counterpart of this model showed no beat phenomenon, indicating that H coverage can improve the robustness of the elementary directions to a certain extent.Figure 8
Effects of H termination (a)~(c), radial (d)~(f), and axial impurity segregation (g, f) on the deflection of the elementary directions. (b) Appearance of the beat phenomenon after removing the H atoms on the surface of a sample in the 7% group (a) and disappearance after rotating the coordinate system by 0.0095π (c). Deteriorated radial symmetry (d) resulting in the beat phenomenon (e) and its disappearance after rotating the coordinate system by 0.0064π (f). Axial symmetry breaking (g) causing no deflection of the elementary directions (h).Although the elementary directions proved robust against the impurity concentration and distribution, different beat periods for the five samples in each group with different concentrations were obtained. This is because once the elementary direction changes, the frequency along this direction changes accordingly. The beat angular frequency is (ωy-ωz). A slight change in ωy or ωz would result in a different beat frequency.The beat frequency can be calculated as follows considering an NW with a fixed cross-section and length [43].
(8)ωy−ωz∝Iy‐Iz,where I is the moment of inertia. However, this equation could not be used to estimate the beat frequency due to the randomness in the distribution of Au atoms. Therefore, symmetry breaking is difficult to quantify. Hence, this study applied an alternative method to characterize the degree of symmetry breaking, i.e., the enveloping cylindrical surface (sharing the same axis with the NW) of the axial mass distribution curve. The radius of the enveloping cylindrical surface in the bare NW (Figure 8(a)) is 1.95 Å and 1.93 Å for the H-covered counterpart. The slightly lower value indicates better symmetry.Then, the effects of the symmetry of the mass distribution were investigated, considering more extreme situations, i.e., when the impurity atoms were radially and axially segregated. Figure8(d) presents a SiNW with Au segregated in half of the region, which was constructed based on the H-covered counterpart in Figure 8(a) obtained by replacing the Au atoms in the other half region with Si atoms. The radius of the enveloping cylindrical surface increased to 2.34 Å, indicating deteriorated symmetry. The excitation direction was the same as above. As observed in Figure 8(e), the model shows a slight beat phenomenon. The auxiliary peak rotates y-z by 0.0064π to obtain the new elementary coordinate system y′-z′. The transformed vibrational signal and spectrum are plotted in Figure 8(f). The rotation angle is also negligible, indicating that radial symmetry breaking slightly affects the deflection of the elementary directions. In comparison, axial symmetry breaking will not cause the deflection of the elementary directions. The model shown in Figure 8(g) is the H-covered counterpart of Figure 9(a) with only one-third of the doped region in the axial direction. The result (Figure 8(h)) shows no beat phenomenon, suggesting that the elementary directions remain unchanged.Figure 9
(a) Quality factor of [1 0 0] and [1 1 0] and (b) beat periodTb SiNWs as a function of concentration.
(a)(b)The above investigations lead to the conclusion that the robustness of the elementary directions in the [1 1 0] SiNWs is hardly affected by impurity Au atoms. However, from a rigorous point of view, the robustness may also be related to the cross-section pattern and boundary conditions (e.g., clamped on one side) or other factors associated with impurity Au atoms, which is a topic for future work. Furthermore, the reason that [1 0 0] SiNWs do not possess this characteristic is also an interesting topic. Actually, for most Au-doped [1 0 0] samples, we failed to find the new elementary axesy′-z′ using Equation (7), irrespective of how we adjust θ, and the FFT curves for y′ and z′ always exhibited an auxiliary peak. Thus, we can preliminarily infer that this behavior is closely related to the symmetry difference and the original elementary directions.This section explains that the unintentional impurity Au atoms will not cause serious deflection of the elementary axes in the [1 1 0] SiNWs. When applied in ultrasensitive resonators, the symmetry breaking induced by the impurity distribution will not produce a large frequency shift. It will not cause disturbance during signal testing. According to the advanced mass detection technique developed by Gil Santos et al. [45], the position and mass of the adsorbate can be obtained by the relative shift of frequencies in two orthogonal directions. Effective restraints on extra frequency shifts caused by the deflection of elementary directions are required. [1 1 0] SiNWs can meet this need because their elementary axes are robust.


## 3.4. Quality FactorQ

Q quantifies the energy dissipation rate of a resonator, defined as the ratio between the total external energy and the average energy loss in one cycle at the resonant frequency. MD can output the total energy, including the external energy (or mechanical energy) and internal energy (or thermal energy). It is inconvenient to extract the external energy to obtain Q; thus, we used an alternative method, i.e., the RMS of the COM [46]. The amplitude of the RMS decreases following the exponential law Ae-ξωt, where ξ is the damping ratio, and ω =2πf0 is the angular frequency, with f0 obtained from FFT. Moreover, Q is equal to 1/2ξ.TheQ values for [1 0 0] and [1 1 0] SiNWs are plotted in Figure 9(a) as a function of concentration. The graph shows an overall downward trend, accompanied by the enlarged error bar. Besides, the impurity distribution considerably magnifies the scattering of Q. The scattering of Q is amplified as the concentration increases. In this case, we cannot safely capture the varying trend of Q as a function of the concentration due to the small sample size. In addition to external damping, atomic friction is also an important factor for energy dissipation [47]. According to Islam et al. [4], impurity atoms can promote friction, which results in decreased Q. Therefore, the downward trend is logically correct.The strong dependence of the beat period,Tb=2π/ωy−ωz, of the [1 1 0] SiNWs on the concentration is shown in Figure 9(b). The longer the beat period, the larger is the deflection of the elementary coordinate system. Although the elementary coordinate system is robust, a slight deflection can lead to an obvious beat phenomenon. From Figures 9(a) and 9(b), one can infer that the deflection of the elementary axes accelerates the vibrational energy dissipation. This feature was not observed in the [1 0 0] SiNWs, probably because of the same geometric symmetry along y and z.More computational resources are required to reveal the relationship betweenQ and the doping concentration than for the other properties, only to obtain the downward trend and increased scattering. Nevertheless, these findings have engineering significance. For example, high Au concentrations should be avoided for standardizing SiNW-based resonators.


## 4. Conclusion

MD was applied to investigate the influences of Au, a commonly used catalytic metal, on the vibrational characteristics of SiNWs, using a newly developed MEAM potential for the Si-H-Au ternary system. The parameterization targets the fundamental properties of H-covered SiNWs with impurity Au, mainly including Au-induced changes in the structure, energy, and elasticity of bulk Si and the properties of H-terminated Si surfaces. Au doping causes evident changes in resonant frequencyf0 and quality factor Q. Both these parameters decrease with the concentration; the former drops linearly, whereas the latter just shows a downward trend accompanied by the amplified error bar. This indicates that Q is sensitive to the distribution of impurities while f0 is not. Besides, the deflection of the elementary axes leads to a higher energy dissipation rate. An interesting observation is that the elementary coordinate system of [1 1 0] SiNWs remains stable. The deflection angle is rather small, even when the impurities are segregated. This observation suggests that when used as ultrasensitive resonators, [1 1 0] SiNWs can effectively restrain the unwanted frequency shift potentially caused by the impurity induced deflection of the elementary axes.



---

*Source: 5790074-2022-09-30.xml*