Multilayer hexagonal silicon forming in slit nanopore

The solidification of two-dimensional liquid silicon confined to a slit nanopore has been studied using molecular dynamics simulations. The results clearly show that the system undergoes an obvious transition from liquid to multilayer hexagonal film with the decrease of temperature, accompanied by dramatic change in potential energy, atomic volume, coordination number and lateral radial distribution function. During the cooling process, some hexagonal islands randomly appear in the liquid first, then grow up to grain nuclei, and finally connect together to form a complete polycrystalline film. Moreover, it is found that the quenching rate and slit size are of vital importance to the freezing structure of silicon film. The results also indicate that the slit nanopore induces the layering of liquid silicon, which further induces the slit size dependent solidification behavior of silicon film with different electrical properties.

In the last decades, low-dimensional materials have gained rapid development along with the microminiaturization trend of semiconductor industry. Graphene, an allotrope of carbon in the form of a two-dimensional (2D), one atom thick, hexagonal lattice, has always grabbed the worldwide attention ever since the experimental preparation in 2004 1 , due to its fascinating properties and nearly infinite applications [2][3][4] . Inspired by the tremendous advancement in graphene, other 2D materials begin to attract an increasingly scientific interest. Recently, strong effort has been invested to provide theoretical and experimental evidence for similar 2D materials composed of other group-IV element, such as silicon 5,6 , which is expected to have a great impact on the development of future nanoelectronic devices. The silicon analogue of graphene was first pointed out by Takeda and Shiraishi and was named silicene by Guzmán-Verri et al. in 2007 7,8 . It has been shown that the silicene with Si atoms packed in a buckled honeycomb lattice, is predicted to possess massless Dirac fermions and to exhibit a detectable quantum spin Hall effect, and other attractive properties 9,10 . Recent experiment obtained an estimated Fermi velocity for the silicene layer of v F = 1.3 × 10 6 ms −1 11 . Due to the specific physical and chemical properties, silicene has broad application prospects in electronic device, energy storage, and so on 12,13 . Currently, Li et al. 14 reported a silicene field-effect transistor operating at room temperature, which opens up new opportunities for 2D silicon for various fundamental science studies and electronic applications.
Different from the carbon, the valence electrons of silicon localized in σ bonds are less mobile than those in graphite, which prevent the formation of extended π electronic states at low binding energies 15 . As a result, the silicon crystal tends to adopt the three-dimensional diamond structure made of sp3-hybridized Si atoms and there seem to be not any natural solid phase of silicon similar to graphite, which precludes the exfoliation method to generate silicene as performed initially in the case of graphene. Now, the most common way to prepare silicene is to deposit silicon on a silver substrate [16][17][18] . Although plenty of theoretical studies had speculated about the existence of silicene, the silicene was not reported experimentally until 2010, when the silicene nanoribbons on Ag(110) with graphene-like electronic signature was observed 19,20 . Thereafter, Daniele et al. 21  tunneling microscopy and scanning tunneling spectroscopy and evidence the presence of corrugated silicene domains.
Though a fair amount of efforts have been devoted to investigating the silicene, the preparation of high quality silicene films is still a major challenge in this field. Recently, some efforts have been devoted to studying the phase transition of confined silicon. Morishita et al. 22,23 used molecular dynamics (MD) simulations to study the formation of nanowire, single-and double-layer silicon in slit pores and the stability is further confirmed by first principles MD calculations up to 300 K, which demonstrate the possibility of the synthesis of novel nanostructures by confinement in nanopores. Afterwards, the electronic and phonon properties of the double-layer silicon computed by Bai et al. 24 suggested that the bilayer hexagonal silicon is a quasi-2D semimetal. In this paper, we perform molecular dynamics (MD) simulations to study the liquid-solid transition of silicon confined in a slit nanopore and provide evidence for the formation of multilayer silicene-like polymorph. The effect of cooling rate and slit size has been considered.

Results
In this study, MD simulations are used to investigate the solidification of liquid silicon confined between two isolated and smooth walls. The liquid silicon is equilibrated at 3000 K first, followed by a cooling process to 300 K, at a quenching rate of 0.1 K/ps. Figure 1(a) shows the potential energy (PE) and volume per atom as a function of temperature. It can be seen that the PE and volume show a strong dependence on the temperature. At the beginning of freezing process, both the PE and the volume drop slowly with the decrease of temperature. At T = 1530 K, the PE suddenly drops by about 0.15 eV/atom from − 3.465 eV/ Atom to − 3.615 eV/Atom. At the same time, the volume per atom increases sharply and the density of silicon film changes from 2.82 g/cm 3 to 2.63 g/cm 3 . To clarify the reason of volume expansion, we plot the size of silicon film in x, y and z directions at different temperature in Fig. 1(b). It can be found that though the thickness of silicon film decreases at T = 1530 K, the plane (x, y) sizes of silicon film increase sharply, which induce the volume expansion of the system.
Usually speaking, such a sudden change in PE and volume suggests the existence of a phase transition [25][26][27] . To better characterize the structural transition during the freezing process, we have plotted the lateral radial distribution functions g xy (r) at different temperature in Fig. 2. It is worth noting that the peak height of the g xy (r) increases significantly with the decreasing temperature, which indicates the structural transition from disordered state to ordered during the cooling process. When the temperature drops to 1500 K, obvious crystal peak begin to appear in the g xy (r), suggesting that the liquid-solid phase transition of silicon film occurs and the system is transformed into ordered crystal structure. With the further decrease of temperature, the crystal peak becomes more obvious and the ordered degree of silicon film increases. In terms of the physical implication of g xy (r), the abscissa value r i of the ith peak always represents the average distance from ith nearest-neighbor shell to the central atom 28 . For example, the r 1 corresponds to the average distance of the first nearest-neighbor shell. Another function of g xy (r) is that scaling behavior of r i /r 1 should be the general feature for the liquid, glassy and crystalline states. At T = 300 K, the abscissa values of the first three peak are r 1 = 2.4 Å, r 2 = 4.15 Å and r 3 = 4.75 Å. The ratios r 2 /r 1 and r 3 /r 1 happen to be about 1.73 and 1.98(approximately 2) respectively, which predicts the solidification structure of silicon film may be composed of hexagonal lattice.
To explore the microscopic details of the structure evolution during the freezing process, simulation snapshots are displayed in Fig. 3. It can be found that the structure has been dramatically transformed during the cooling process. As shown in Fig. 3(a-c), random silicon island appears in the liquid silicon before the liquid-solid phase transition. However, it is worth noting that the silicon islands in this stage can persist only for a short time and are easy to break into disordered structure inversely. If the size of islands exceeds a critical value, the islands can stable exist as a nucleus as shown in Fig. 3(d), which indicates the beginning of liquid-solid phase transition. With the further decrease of temperature, the liquid silicon is simultaneously nucleated at individual islands. Subsequently, the islands begin to grow outwards by incorporating the disordered structure. Then, the carbon islands will encounter each other and combine to form a complete polycrystalline film. Because the lattice orientation of every island is quite different, there may be some grain boundaries in the solidification structure. Eventually, the 2D silicon exhibits a bilayer hexagonal polycrystalline structure as shown in Fig. 3(i), agrees well with previous results 24 .
The phase transition can also be indicated by a sudden change in coordination number. In this study, the coordination of each atom counts neighbors within a radius of 3 Å. As shown in Fig. 4, it can be found that the coordination fractions are strongly dependent on the slit size. Before solidification, the average coordination first increases and then decrease slowly with the decreasing temperature. In the stage of liquid-solid phase transition, the average coordination varies dramatically from 4.45 to 4.02. At this point, the fraction of 4-fold coordination increases sharply while the fraction of the other coordination types decreases steeply. At T = 300 K, the fraction of 4-fold coordination reaches nearly 100 percent, rather different from the 3-fold coordination dominated graphite, although the silicon film and graphite both have hexagonal lattice. This is because the silicon atoms between two layers are combined by σ bonds rather than the weaker π bonds 15 .
Next, we will discuss the freezing process of silicon film at three different quenching rates (QRs). As shown in Fig. 5(a), the QR has little effect on the change of PE before solidification. With the decreasing temperature, there is a sharp decline in PE curve, indicating the liquid-solid phase transition. It is worth noting that the decline ranges of PE curves decrease and are getting less stark as the QR increases. At QR = 0.1 K/ps, the PE suddenly drops at T = 1530 K, indicating the completion of solidification at this temperature. However, the PE curve does not show a dramatic decline in a certain temperature but in a range of temperature when the QRs are 1 K/ps and 10 K/ps, usually indicating the generation of amorphous structure. To further discriminate the solidification structure, the g xy (r) of solidification structure at different QRs are plotted in Fig. 5(b). It can be found that the crystal peak of g xy (r) at QR = 0.1 K/ps is very obvious, suggesting that the system is dominated by ordered structure. At QR = 10 K/ps, the second peak of g xy (r) splits into two subpeaks and there is no crystal peaks, which indicates that the solidification structure is amorphous. At QR = 1 K/ps, the crystal peak and the splitting second peak exits in the g xy (r) at the same time, indicating the coexistence of ordered and disordered structures.
To explore the microscopic details of the structure evolution at different QRs, Fig. 6 shows the simulation snapshots at different QRs. As shown in Fig. 6(a), the nucleation and growth of liquid silicon begins with four stable islands at QR = 1 K/ps. However, due to the fast cooling process, the crystal nuclei do not have enough time to grow up and a part of liquid silicon directly freezes into amorphous structure before converting to ordered structure, which induces the coexistence of ordered and disordered structures. This phenomenon becomes even more clear at QR = 10 K/ps. As shown in Fig. 6(d-f), there are only a few silicon around the nucleation island transformed into ordered structure and most of liquid silicon directly freezes into amorphous structure. As a result, the g xy (r) shows the amorphous characteristics.
As clarified above, the QR is of vital importance to the freezing structure of silicon film. Next, we will investigate how the slit size affects the solidification behavior of the liquid silicon at QR = 0.1 K/ps.   Figure 7 shows the PE per atom as a function of temperature for different silt sizes D. It is can be found that there is obvious sudden drop in the PE curves at T = 1530 K and 1360 K for D = 11 Å and 13.2 Å respectively, suggesting that the silicon films solidify at that temperature and is dominated by the ordered structure. It is also worth noting that the freezing point of silicon shifts towards the lower temperature when D changes from 11 Å to 13.2 Å. For D = 12.3 Å and 15.9 Å, the PE curve does not show a dramatic decline in a certain temperature but in a range of temperature, indicating the generation of amorphous structure.
To further study the microscopic structure of silicon film after solidification, simulation snapshots at 300 K are investigated, as displayed in Fig. 8. It can be found that the solidification structure is very different for different slit sizes. The silicon film for D = 11 Å and 13.2 Å is composed of two and three atomic layers with hexagonal lattice, respectively. At D = 11 Å, the atoms in the two layers are mainly 4-fold coordinated. At D = 13.2 Å, the atoms in the outer layers are mainly 4-fold coordinated, while the atoms in the inner layers are mainly 5-fold coordinated, so that the whole film is a 4-5-4 coordinated sandwich structure. At D = 15.9 Å, the silicon films is composed of both ordered hexagonal structure and disordered structure. Here, the hexagonal structure is constructed by four atom layers coordinated as  4-5-5-4. Different form the case at D = 15.9 Å, though the silicon films also contains both ordered hexagonal structure and disordered structure at D = 12.3 Å, the hexagonal structure is not plane but buckled, as shown in the inset of Fig. 8(d). The buckled hexagonal structure is induced by the protuberance of one 3-fold coordinated atom in every ring. The protruding atom lines arrange alternatively in the two layers, which induces the fluctuation of silicon film.
In order to clarify the reasons for such the slit size dependent solidification behavior of silicon film, we have plotted the density distribution functions along the confined direction for the liquid (1800 K) and solid (300 K) silicon in Fig. 9. The results clearly suggest a structural correlation between the liquid and solid silicon. At D = 11 Å, 13.2 Å and 15.9 Å, the density distribution functions of the liquid and solid silicon both have two, three and four peaks as shown Fig. 9, indicating the two phases are both composed of two, three and four atomic layers, although the density peaks of liquid silicon are rather smooth than that of the solid one. At D = 12.3 Å, each main density peak of the liquid silicon split into two subpeaks,  which becomes more obvious in the solid silicon corresponding to the buckled hexagonal structure. Based on the above, one can draw a conclusion that the structural difference of silicon film has already appeared before solidification and the layer number is determined by the slit size. With the decreasing temperature, the layering phenomenon is increasingly evident. This result is in good agreement with previous experimental and numerical studies that confined liquid undergoes a layering transition nearing the solid wall [29][30][31] . Therefore, it can be deduce that the slit nanopore induces the layering of liquid silicon, which further induces the slit size dependent solidification behavior of silicon film.
The stability of stand-alone multi-layer silicon (without the confinement) has been examined by first-principles MD calculations. The temperature and pressure of 2-layer (64 atoms), 3-layer (96 atoms) and 4-layer (128 atoms) structures are set to be 400 K and zero respectively, using the Nosé-Andersen method. The time step is 1fs. Figure 10 shows the snapshots of the first-principles MD simulation at 10 ps. It is found that there is no obvious structural transition except for slight vibrations of silicon atoms, which indicates that the stand-alone bilayer silicon can still stably persist at 400 K. After confirming the stability of the multi-layer hexagonal silicon polymorph, the electronic properties are also explored by first-principles MD calculations. The Brillouin zone was sampled with (4 × 4 × 1) k points of a Monkhorst-Pack grid. Figure 11 shows the band structures of multi-layer hexagonal silicon after full geometry optimization. It can be found that the 2-layer silicon shows a classic semimetal characteristic. With the increase of layer number, the silicon tends to become a metal, as shown in Fig. 11(c). Figure 11(d) shows that the buckled 2-layer silicon is semi-conductor of electricity.

Discussion
MD simulations have been performed to study the cooling process of liquid silicon confined in slit nanopore. With the decrease of temperature, the silicon island randomly appears in the liquid first and then grows up to grain nucleus, indicating the start of freezing transition. Subsequently, the liquid silicon transforms to multilayered hexagonal film in a short period, along with dramatic change in potential energy, atomic volume, coordination number and lateral radial distribution function. For the ordered hexagonal structure, the atoms in the layer adjacent to the wall are mainly 4-fold coordinated while inner layer are dominated by 5-fold coordinated atoms, which is rather different from the multilayer graphene.
The quenching rate and slit size are both very important to the solidification behavior of silicon film. On the one hand, the increasing quenching rate will lead to the non-crystallizing of the liquid silicon. On the other hand, the increasing slit size not only affects the structure of each layer, but also introduces new layers in the multilayer film. Our results clearly show that the slit nanopore induces the layering of liquid silicon, which further induces the slit size dependent solidification behavior of silicon film. First-principles MD calculations are used to confirm the structural stability of multi-layer silicon. The band energy calculations suggest that layered silicon with different structures has different electrical  properties. These findings provide physical and dynamic insights into the liquid-solid phase transition in 2D silicon and predict a possible way to fabricate the 2D silicon materials with different properties.

Methods
In this study, MD simulations are performed to investigate the solidification of liquid silicon confined between two isolated and smooth walls. The Stillinger-Weber empirical potential is employed to describe the silicon-silicon interaction 32 , which comprises two and three body terms, energetically favoring tetrahedral coordination of the atoms. We utilize the 12-6 Lennard-Jones (LJ) potential, with a well depth ε = 0.310 kcal/mol and size parameter σ = 3.804 Å [33][34][35] , to describe the silicon-wall interaction. Although the SW potential has known limitations 36 , it provides a reasonable description of the crystalline and liquid phase [37][38][39] and has been widely used to study the structure transition of silicon [40][41][42] .
For all simulations, we adopt the constant number, lateral-pressure, temperature (NPT) ensemble and use the velocity-verlet algorithm with an integration time step of 1 fs. In this work, we use the MD package LAMMPS to perform the simulations 43 . The periodic boundary conditions are applied only in the lateral directions parallel to the wall at zero pressure, while in the vertical direction perpendicular to the wall, the slit size D between two walls is originally set to be 11 Å. For each D, the liquid silicon is equilibrated at 3000 K first, followed by a cooling process to 300 K, at a quenching rate of 0.1 K/ps. Then, the final structures are obtained after a relaxation process for 200 ps at 300 K. To confirm the structural stability and study the electrical property of multi-layer silicon, first-principles MD calculations are used based on density functional theory within the generalized-gradient approximation (GGA) 44 with the Perdew-Burke-Ernzerhof (PBE) functional and ultrasoft pseudopotential 45 .