Molecular dynamics simulation of pressure effect on silicene nanoribbons

Introduction: In 2D materials, the boundary of silicene formed as nanoribbons plays an essential role in synthesis and can be controlled to achieve different characteristics. Our study aims to investigate the structural preference depending on the pressure tolerance and boundary dependency. Methods: The main methodology used in our study is molecular dynamics simulation with Stillinger Weber potentials. Our simulation was carried out on 2D models of honeycomb silicene obtained through high and low pressurized cooling from the liquid state and then heat annealing for a decent time. The final configuration of silicene will be investigated in terms of structures and thermodynamic properties. Results: We found that the effect of the cooling process under high pressure formed a 4-fold ring structure, while at low pressure, 2D honeycomb networks were recovered but with different degrees of defects depending on the boundary condition. The main difference between several transitions of 2D silicene is discussed via the evolution of total energy and the change in coordination number and bond-ring distribution. Conclusions: This study provides insights into the dependency of the structure of silicene on the pressure and boundary, repre-sentedby the first-order transition at lowpressure and a congregationof disorderedlow-numbered rings into the ordered tetragonal formation at high pressure. Notably, our results have shown that silicene nanoribbon materials can be controlled by pressure to obtain unfamiliar structures such as pentagonal and tetragonal networks.


INTRODUCTION
For material science in the last decade, twodimensional groups of materials led by graphene have attracted extensive attention from our scientific community. The properties of graphene have promoted many new studies in 2D materials in which silicene has become a potential candidate for developing semiconductor generation 1,2 . Unlike the flat monolayer in graphene, silicene has a buckling structure. In addition, silicene preserves all the advantages of bulk silicon materials as its electronic configuration, and isotropic affinity can easily be put into practice due to compatibility with current semiconductor technology. In addition, silicene has several outstanding properties in applications such as ferromagnetic, semimetallic, quantum Hall effect, giant reluctance, optoelectronics, spin electrons, and superconducting materials [1][2][3][4][5][6] . Many 2D structures of silicene and Si configurations, such as SiC, silicene bilayer, and defective silicene, have been studied recently [7][8][9][10] . The free-standing 2D structure of silicene associated with amorphous and rapid cooling liquid has been studied in theory and experiment [11][12][13] . Amorphous silicene was found to form through the fast quenching rate in Ref. 13. The results have shown that the 2D amorphous structure of silicene is very complex, including many types of rings, mainly from 3-fold rings to 6-fold rings. In addition, a study by Deb et al. showed a phase transition of amorphous silicene due to the pressure effect between the low-density and high-density phases accompanying the liquid-amorphous phase transition from approximately 900 K to 1300 K 14 .
Recently, H.A. Huy et al. studied the phase transition from the disordered amorphous/liquid state of 2D silicene to the "Cairo-tiling" pattern of penta-silicene and the 4-fold ring structure of tetra-silicene at high density by the MD simulation method 15 . These results estimate the critical pressure while the density of phase changes and can be used as a reference for future silicene development. This conclusion was also confirmed by M. Qiao et al. in 2017 16 . In 2019, V.V. Hoang et al. showed that due to very weak interactions with the walls, the structural features of the final configuration of tetra-silicene should be identical to those of free-standing tetra-silicene 17 . In addition to hexa-silicene materials, the cooling process and thermal stability of penta-silicene materials have also been studied under the high-pressure effect of the amorphous state 18 . The thermal stability of the models demonstrates that materials can be controlled by pressure to obtain the desired structure. In experiments, silicene is synthesized via epitaxial growth as nanoribbons, which are called SiNRs 19 . In nanoribbon formation, the size effect has been considered an important factor in controlling many characteristics. Due to the confinement of the charge carrier in the lowdimensional structure, manipulating the ring size and thus the band gap could be feasible. In fact, the electronic transport study of graphene nanoribbons by M.Y. Han et al. revealed that the energy gap scales inversely with the ribbon width 20 . Recently, an experiment discovered that single-and double-strand SiNR formation could comprise pentagonal rings 21 . This is similar to the penta-graphene formation growth on Cu substrates found by K. Xia et al. in 2016 22 . The different formations of the 2D structure could involve a variation in the strain applied (which is considered equal to the pressure of the model). Therefore, our study aims to investigate the variation in SiNR structure under low-and high-pressure conditions.

MATERIALS AND METHODS
This study is carried out via the molecular dynamics (MD) method to model SiNRs. The potential of Si-Si interactions was proposed by F. H. Stillinger & T.A. Weber 23 , where the pair and three-atom interaction is expressed as: The first two-body term U 2 consists of ionic repulsion and attraction for short-range interaction: The second term U 3 is a summary of all triplets of a three-body interaction: Details of the SW potential, including the parameters, can be found in Ref. 23 . SW potential has been widely used for the simulation of Si or Si-based systems for the last 30 years. In the study of 2D silicene, the SW potential has several advantages, such as being parameterized easily and appropriately determining the melting/crystallization temperature. Recent studies of our simulated silicene cooling models from the melting state using SW potential give good results 13,15,18,24,25 . Models of 6-fold ring structures with both zigzag and armchair boundaries are investigated in two cases: low-pressure (0.31 Mbar) and high-pressure (0.73 Mbar). The simulated models consist of 10 4 atoms, and the details of the relaxed size for the four models are presented in Table 1. The initial temperatures of the models are 50 K. As shown in Figure 1, these models are melted isovolumetric with melting rate γ = 2 × 10 11 K/s until the temperature of the models reaches 3500 K and then relaxed for a while until they completely melt. After that, the models are cooled with the same rate γ = 2 × 10 11 K/s down to 300 K and then relaxed in a time t = 10 ns corresponding to the heat annealing process to investigate the stable province over time. The simulation process is compared under two different pressure conditions to determine the effect of pressure on the formation of the material structure. The MD simulation process was calculated via the LAMMPS package based on the Nose-Hoover barostat implemented for controlling the pressure 26 . The ensemble conditions NVT are applied for lowpressure models, while NPT conditions are used in high-pressure simulations. For nanoribbon modeling, periodic boundary conditions (PBCs) are applied only in the x Cartesian directions. In the z and y directions, a fixed boundary with elastic reflection behavior is applied. The temperature is corrected via simple velocity rescaling, and the Verlet algorithm is employed. For analysis, we employed ISAACS software with the 'shortest path' rule to calculate the ring statistics 27 . The selected cutting radius is r cuto f f = 2.85 Å in calculating the coordination number, bond-angle distribution, and interatomic distance distributions of the models. The cutoff radius is determined by the position of the first minimum after the first peak in the radial distribution function (RDF) of models obtained at the final compression process. VMD software is used for 2D visualization of atomic configurations 28 . The results are averaged over two independent runs to improve the statistics.

Results of the SiNRs after solidification under different pressure conditions
First, we present the final result of SiNRs manifested from the solidification process under two different pressure conditions. Figure 2 describes the visual images of the structure obtained after cooling and relaxing at 300 K of SiNRs in two boundaries under different pressure conditions. At low pressure, both the armchair and zigzag SiNRs tend to recover the honeycomb structure in Figure 2a and c; however, the armchair model exhibits a less homogeneous ribbon structure with several domains of the pentagonal ring. At high pressure, the honeycomb structure will be replaced by the 4-fold ring formation or the pentagonal ring, as shown in Figure 2c, d-1, and d-2. From the ring formation in Figure 2d-1, the zigzag model could manifest several domains of the "Cairo-tiling" pentagonal ring sandwiched between two 4-fold ring formations. Another possible location for the "Cairotiling" pentagonal ring is the grain boundary, as seen in Figure 2d-2. Hence, we show that when cooled under high pressure, the armchair model appears to have a more homogeneous 4-fold ring structure than the zigzag model. Insights into the dependency of SiNRs on the width boundary are discussed in the next section.

Detailed analysis of the SiNR structure
To study the structural details of the models, we investigate the radial distribution function, interatomic distance, and bond-angle distribution of all models at 300 K. The RDF describes the distribution of bonding distances between Si atoms from 0 Å ÷ 10 Å. Figure 3a shows the RDF of the models at 300 K when the structural models stabilize. It shows a significant difference between the models. In the high-pressure models, the distances of the first, second, and third peaks are r 1 = 2.44 Å, r 2 = 1.4r 1 , and r 3 = 2r 1 , while in the lowpressure models, they are r 1 = 2.33 Å, r 2 = 1.73r 1 , and r 3 = 2r 1 . From Table 2, the high-pressure models correspond to a 4-fold ring structure, while the lowpressure SiNRs still maintain the honeycomb structure of 6-fold rings. The average Si-Si bonding characteristics in Figure 2b show that the impact of pressure on the structural ring formation in the nanoribbon is noticeable. Similarly, the peaks of the Si-Si interatomic distance in the high-pressure models shift to a higher value, indicating sp 2 ↔ sp 3 hybridization from 3-coordinated to 4-coordinated rings. Meanwhile, the interatomic distances of low-pressure SiNRs are slightly higher than those of pristine silicene, which has a Si-Si bond at 2.28 Å. We also indicate that the boundary affects the distribution of interatomic distances in SiNRs at 300 K. The broadening of interatomic distances in low-pressure zigzag SiNRs shifts at higher values, expressing less defective 6-fold ring formation. In high-pressure models, zigzag SiNRs have a broader distribution, which means that their ring structure is less uniform than that of armchair SiNRs. Figure 4 shows the bond-angle distribution of the final models. In the low-pressure models, the distribution with a sharp peak at 118 • is close to the value of 116.3 o of the 6-fold ring ideal models. However, this distribution is quite broad, ranging from 90 • to 151 • in the zigzag boundary model, while the armchair boundary model has a wide distribution from 47 • to 180 • . The distribution ratio at the sharp peak of the zigzag boundary model is higher than that in the armchair boundary model. A broad distribution represents structural defects in the models, so the zigzag boundary model has fewer defects than the armchair model. Similarly, the bond-angle distribution in the high-pressure models has two sharp peaks at 88 • and 162 o . This value features the 4-fold ring structure. In addition, the graph of the bond-angle distribution in the high-pressure models also shows that the armchair boundary model is less defective (ratio of higher peaks and fewer minor peaks). The broad angular distribution is associated with the defective pentagonal structure of "Cairo-tiling" in both low-pressure   and high-pressure models. Our findings correspond to the confirmation of the pressure effect on the structural model in our previous discussion. We also note that this pressure effect is similar to what is found in the compression of 2D amorphous silicene 18 and confined silicene from the melt 17 . Furthermore, the boundary dependency of the SiNR structure is considerable and could lead to different arrangements at low and high pressures. The obtained structures from Figure 2 are consistent with the RDF function and bond-angle distribution results in Figure 3 and Figure 4. From the overall features of the RDF and bonding analysis above, the SiNR structure can be rearranged into several morphology patterns depending on the boundary and high-pressure conditions. Indeed, our simulated SiNRs with a 6-fold ring structure recovered at low pressure with different defective degrees. Meanwhile, the 4-fold ring structure is formed at high pressure. Interestingly, a significant amount of "Cairotiling" pentagonal structure could manifest in the zigzag models in two models.

Thermodynamic properties of SiNRs in the cooling process
To obtain more insight into the phase transition, we calculated the total energy of the systems during the cooling process. The dependence of the system's total energy on the temperature during the cooling process at low-pressure and high-pressure conditions is shown in Figure 5. We found that at 300 K, the lowpressure model has lower total energy than the highpressure model in both armchair and zigzag boundaries. This tendency indicates that the model cooled at low pressure undergoes a structural transition, while the high-pressure model does not express any significant change. The curves of the total energy of both armchair and zigzag boundaries at high pressure are less disrupted compared to the low-pressure transitions. Here, we also found that the different boundaries in SiNRs could lead to two distinct cooling processes. As shown in Figure 5a and Figure 5b, there is a remarkable difference in the energy variation and the phase transition temperature between the armchair and zigzag boundaries. In the armchair boundary, the phase transition temperature T C = 1400 K, while in the zigzag boundary, T C = 1550 K with the lowpressure condition and T C = 1500 K with the highpressure condition. Regarding the structures in Figure 2 that obtained the same cooling process from the melting state, the energy evolution of the four SiNR models will indicate the general dependency of the applied pressure on the final preferred formation manifested at low temperature. The energy curves suggest that the low-pressure structural transition is a firstorder transition driven by the transformation of disordered low-numbered rings into a stable honeycomb structure. Meanwhile, the high-pressure model is a congregation of disordered low-numbered rings into a fine-ordered tetragonal formation. More detail on the SiNRs structure during their solidifications will be analyzed via the coordination number and bond-ring distribution of models in the cooling process, as shown in Figure 6 and Figure 7. Under low-pressure conditions, the 3-coordinated ratio of the armchair model is approximately 80%, and in the zigzag model, it is 90%, with some ratios of 4coordinated and 2-coordinated. Similarly, most of the atoms participate in the 6-coordinated structure (Figure 7a, Figure 7c) with a ratio of over 70% in the zigzag model and approximately 45% in the armchair model. In the models that are cooled under low pressure, both  boundary forms are obtained in 3-coordinated models with a high ratio. However, the ring distribution expresses a significant proportion of the 5-fold rings (40%) in the low-pressure armchair boundary model, indicating that this structure has a complex formation. Meanwhile, the low-pressure zigzag boundary model only contains approximately 21% of the 5-fold ring. These 5-fold rings take the shape of "Cairotiling" first introduced in 2D ice sheets by J. Chen et al. in 2016 29 and are similar to the structure of compressed amorphous silicene by H. A. Huy et al. in 2019 18 . As seen from the result in Figure 1a, there is a domain of "Cairo-tiling" pentagonal ring formation sandwiched between two honeycomb domains. In Figure 2c, the presence of pentagonal rings only exists as defects in the low-pressure zigzag boundary model. Therefore, the low-pressure model obtained from the zigzag boundary is more homogenous than the armchair boundary model. In high-pressure conditions, the coordination numbers in the zigzag and armchair boundaries also differ. The 4-coordinated ratio of the armchair model at 300 K is nearly 90%, while the ratio of the zigzag model is approximately 70% and appears 2-coordinated in this model. This difference indicated that the model cooled under high pressure will obtain high 4-coordinated ratios, which is characteristic of the 4-fold ring structure. Indeed, Figure 7b and Figure 7d show that the number of atoms participating in the 4-fold ring structure in the armchair model is 100%, and in the zigzag model, it is 75%. Therefore, the high-pressure model is formed as a gradual increase of the 4-fold rings in the liquid state (composed as a mixture of 3, 4, and 5-fold rings) into the dominant 4-fold ring formation. In contrast, during cooling at low pressure, the 3-and 4-fold rings are diminished and replaced by the 5-and 6-fold rings, resulting in two structures, as shown in Figure 2a and c. Hence, the pressure condition triggers sp 2 ↔ sp 3 in 2D silicene, similar to what occurs in silicon 14 and amorphous silicene 15 . Each condition favors the intermediate pentagonal ring as defective domains or grain boundaries depending on the initial ribbon boundary. This phenomenon provides a hint into the fabrication of unfamiliar 2D structures by controlling the pressure conditions and emphasizes the role of the ribbon boundary in the detailed structure of 2D silicene.

CONCLUSIONS
We studied the material forming process at different pressures by the MD simulation method with SW interaction potential. We found that the armchair and zigzag boundaries have distinct effects on the final configuration. The calculated results also show that the 6-fold ring structure still contains different degrees of defective domains, especially if the initial boundary is armchair. At high-pressure cooling, the 4-fold ring structure is approximately 0.74 Mbar, with the armchair boundary model exhibiting an almost ideal 4-fold ring structure. The thermodynamic evolution and ring distribution analysis imply that low pressure is a first-order transition from low-numbered rings to attain the well-known honeycomb network. However, at high pressure, the lownumbered rings are likely to aggregate into a fine tetragonal structure. Therefore, the results of the preference of the SiNR structure for boundary and pressure conditions could have directional meaning in fabrication. Controlling the material structure under the effect of pressure could be useful for the practical synthesis of different potential nanoribbon structures.