Defects in Silicene: Vacancy Clusters, Extended Line Defects, and Di-adatoms

Defects are almost inevitable during the fabrication process, and their existence strongly affects thermodynamic and (opto)electronic properties of two-dimensional materials. Very recent experiments have provided clear evidence for the presence of larger multi-vacancies in silicene, but their structure, stability, and formation mechanism remain largely unexplored. Here, we present a detailed theoretical study of silicene monolayer containing three types of defects: vacancy clusters, extended line defects (ELDs), and di-adatoms. First-principles calculations, along with ab initio molecular dynamics simulations, revealed the coalescence tendency of small defects and formation of highly stable vacancy clusters. The 5|8|5 ELD – the most favorable extended defect in both graphene and silicene sheets – is found to be easier to form in the latter case due to the mixed sp2/sp3 hybridization of silicon. In addition, hybrid functional calculations that contain part of the Hatree-Fock exchange energy demonstrated that the introduction of single and double silicon adatoms significantly enhances the stability of the system, and provides an effective approach on tuning the magnetic moment and band gap of silicene.

Defects are almost inevitable during the fabrication process, and their existence strongly affects thermodynamic and (opto)electronic properties of two-dimensional materials. Very recent experiments have provided clear evidence for the presence of larger multi-vacancies in silicene, but their structure, stability, and formation mechanism remain largely unexplored. Here, we present a detailed theoretical study of silicene monolayer containing three types of defects: vacancy clusters, extended line defects (ELDs), and di-adatoms. First-principles calculations, along with ab initio molecular dynamics simulations, revealed the coalescence tendency of small defects and formation of highly stable vacancy clusters. The 5j8j5 ELD -the most favorable extended defect in both graphene and silicene sheets -is found to be easier to form in the latter case due to the mixed sp 2 /sp 3 hybridization of silicon. In addition, hybrid functional calculations that contain part of the Hatree-Fock exchange energy demonstrated that the introduction of single and double silicon adatoms significantly enhances the stability of the system, and provides an effective approach on tuning the magnetic moment and band gap of silicene. S ilicene, a hexagonal mesh of silicon atoms, has attracted significant interest from both academia and electronics industry because it exhibits comparable electronic characteristics as graphene, e.g., the existence of the Dirac cone and the realization of the quantum spin Hall effect [1][2][3][4][5] . Besides that, silicene could be synthesized and processed using mature semiconductor techniques, being more readily integrated into existing electronics than graphene. Unlike graphene, silicene forms a slightly buckled monolayer structure 6-7 that allows more flexibility to tailor their band structures and functions [8][9][10][11][12] . Despite such single-layer silicon has been theoretically speculated for two decades 6 , the actual fabrication of silicene on many substrates, such as Ag(111) [13][14] , ZrB 2 (0001) 15 , (2 3 1)-reconstructed Au(110) 16 , and Ir(111) surfaces 17 , has been reported only very recently. The successful synthesis of silicene ribbons and sheets shows a promising application potential in nanoscale materials and devices 18 .
Defects are often the first concern in the real application of monolayer materials. Vacancy defects, which are readily induced by laser irradiation and electron beam, are almost inevitable in the fabrication and processing of monolayers, and sometimes, small defects are introduced purposively for specific applications [19][20] . For example, the creation and elimination of point defects have been shown to provide a powerful avenue to tune the local structure, thermal stability, and band gap of low-dimension materials [21][22][23][24] . In addition, the two dimensional materials with defects has been proposed to be a superior membrane for gas separation [25][26] .
Several theoretical works have concentrated on the structure and energetics of point defects in free-standing silicene. For example, Ö zçelik et al. 23 studied the atomistic mechanisms in the self-healing of vacancy defects, and uncovered that the reconstruction process is driven by the energy gained through new Si-Si bond formation. Guo et al. 22 systemically investigated the structures, formation energy, and migration behaviors of typical single vacancy and double vacancies. Berdiyorov and Peeters 27 reported that vacancy defects include local structural changes and considerably reduce the thermal stability of silicene. These prior studies, generally focusing on the single, double, and triple defects in silicene, deepen one's understanding on the formation and diffusion mechanisms of vacancies in silicene. However, very recent experimental studies have provided strong evidence that a large number of defect clusters are present in the scanning tunneling microscope images, with more than three missing silicon atoms 28 . Although the presence of multi-atom vacancies is of fundamental importance for understanding the formation and operation of layered materials, their configurations, stabilities, and formation mechanisms are far from being fully comprehended.
In this contribution, we performed density-functional theory (DFT) calculations to investigate the reconstruction, coalescence, and diffusion behaviors of vacancy defects in free-standing silicene. The structure and energetic of vacancy clusters formed by removing 1-6 atoms were systematically considered, where upon coalescence occurs spontaneously to reduce dangling bonds on the edges. We then focused on five architectures of extended line defects (ELDs) embedded in silicene, showing that ELDs have lower formation energies compared to vacancy clusters, which arises from the reconstruction of periodic defects. Finally, we studied the adsorption and diffusion of silicon atoms on silicene layers, and uncovered that such adsorption process can enhance the stability of the system and open sizeable band gaps in silicene.

Results
The formation energy per unit-cell involved in the reconstruction, E f , which indicates the stability of a defect in silicene, can be determined as where E T is the total energy of a silicene plane containing defects. N denotes the number of atoms in the unit cell, and E gr is the energy per silicon atom in a pristine silicene sheet. Note that the magnitude of E f depends on the number of missing atoms, g, in the silicene sheet. In order to compare the relative stability of systems containing different numbers of vacancies, we further compute the formation energy per missing atom, E f 9, which is defined as Reconstruction and coalescence of vacancies in silicene. We started by considering all possible configurations for the smaller vacancy clusters (1-3 missing atoms), and the most-likely configurations for the larger clusters (4-6 missing atoms). These vacancy defects were created through removing 1-6 silicon atoms from a pristine lattice. We noticed that even a single vacancy defect could lead to a significant local distortion of silicon's hexagonal arrangement (see Fig. 1a), in agreement with previous literature 22 .
Figures 1e-1j illustrate the initial structures of silicene monolayer with 4-6 vacancies and their fully reconstructed counterparts. Apparently, more complicated defect configurations are formed upon the removal of more than three silicon atoms from the silicene sheet. The V 6 (18 6 ) shown in Fig. 1i is found to be the most unfavorable structure, which is presumably because of the presence of the under-coordinated Si atoms. The same observation, i.e., undercoordinated Si atoms destabilizes systems, is also applied for the V 3 (5-10 1 -5), V 4 (5-12 2 -5), V 5 (5-16 5 ), and V 5 (5-11 1 -5) vacancies. The number of dangling bonds largely determines the relative stability of the layer with the same missing atoms. For example, the V 4 (5-5-5-9) structure is more stable than V 4 (5-12 2 -5), due to existence of two dangling bonds (the subscript in parenthesis) in www.nature.com/scientificreports SCIENTIFIC REPORTS | 5 : 7881 | DOI: 10.1038/srep07881 the latter case. Our simulations also demonstrated that instantly removing a number of atoms from silicene causes the bending or rippling of the monolayer so as to considerably reduce its surface area. Indeed, the rippled structure decreases the energy of the systems and may stabilize vacancies in the silicene layer.
To assess whether the proposed vacancy defects in silicene is locally stable, we continued to perform ab initio molecular dynamics (AIMD) simulations on the reconstructed structures shown in Fig. 1. During a 5 ps AIMD run at 500 K, we observed that the metastable V 4 (5-12 2 -5), V 5 (5-16 5 ), and V 6 (18 6 ) structures destabilized and reconstructed again (''re-reconstructed'') to reduce all dangling bonds in the system (see Fig. 2). Importantly, these ''rereconstructed'' structures remain less stable than the V 4 (5-5-9-9), V 5 (5-11 1 -5), and V 6 (5-5-10-5-5) shown in Fig. 1. For example, the computed formation energy of V 4 (5-5-5-9) is 0.09 and 0.70 eV lower than that of r-V 4 and V 4 (5-12 2 -5), respectively, confirming the validity of our DFT-predicted vacancy structures. We now discuss the coalescence of small vacancies in silicene. The tendency of coalesce small defects into larger ones can actually be roughly estimated by the analysis of the E f 9 values shown in Fig. 1. For example, the formation energy of two V 1 (5-5) is significantly larger than that of one V 2 (5-8-5), suggesting that divacancy is more preferable compared to two single vacancies in the silicene. Our AIMD runs at 500 K further confirmed the above hypothesis: two adjacent single vacancies indeed coalesce and transform to the V 2 (5-8-5) structure at 2 ps (see Fig. 2d).
The above conclusion, i.e., small defects tend to coalesce into larger ones, can also be achieved for five-atom and six-atom vacancies. Figure 3 shows the six possible initial structures in the relatively larger (8 3 8) Fig. 3c]. After optimization we found that V 5 has the lowest energy whilst V 2 1 V 3 is the most unstable structure, suggesting that the two separate defects are inclined to merge into a larger vacancy cluster. Similar findings can also be achieved in the case of 6-atom vacancies. As shown in Fig. 3, the three separated V 2 defects have the highest energy, and when the three V 2 defects are arranged in a line and adjacent, the energy decreased and then the three V 2 defects coalesce into one six-atom cluster vacancies. The coalescence process in silicene is probably driven by the weak but long-ranged interactions between the vacancies, due to electronic effects 29 . Since long-range correlations are expected to be more important in anisotropic and low-dimensional structures 30-32 , we recalculated the E f values for the V 1 (5-5), V 2 (5-8-5), and V 6 (5-5-10) structures using the dispersion-inclusive DFT-D2 method 33 . For small vacancy clusters like V 1 , the E f values are almost identical no matter whether the van der Waals (vdW) interactions are taken into account in our DFT calculations. Notably, vdW forces contribute more to the bonding in larger vacancy clusters. As   exemplified in the V 6 (5-5-10) structure, the computed formation energy from DFT-D2 is 8.4% larger than that from the PBE functional. However, our systematic tests have confirmed that the inclusion of the vdW energy does not change the relative stability of different defect systems.
Structures and stability of extended line defect in silicene. ELDs in a two-dimensional (2D) lattice can be seen as a relative displacement of the two halves of the lattice, with possible removal of certain atoms located on or adjacent to the dislocation. In the case of silicene, the 5j7 and 4j8 ELDs are formed along the zigzag and armchair boundary, respectively (see Figs. 4a and 4b). After full relaxation, the neighboring half-lattices couple to each other and form periodic arrays of vacancies. The 5j8j5 ELD shown in Fig. 4c contains two pentagons spatially separated by an octagonal ring, which is obtainable by removing silicon dimmer stilted 30u from the zigzag chain.
In analog to graphene [34][35][36] , array of defects can be reconstructed from divacancies plus a Stone-Wales (SW) transformation. This leads to the formation of a double-pentagon double-heptagon structure (d5d7; Fig. 4d), and a triple-pentagon triple-heptagon structure (t5t7; Fig. 4e). Note that the above two ELDs, which can be obtained by irradiation, have been demonstrated to be stable topology for the reconstruction of an isolated divacancy in graphene 36 . Given that the buckled silicene requires a lower SW transformation barrier than graphene 24 , one would reasonably expect the existence of the d5d7 and t5t7 ELDs in the silicene monolayer.
The relative stability of ELDs in silicene can be estimated by the formation energy per unit, E f 9, which is defined as where d is length of the supercell along the ELD direction. Our calculations show that the 5j8j5 defect is the most stable structure, with a formation energy 0.06 eV/Å lower than that of the t5t7 defect. We have also compared the formation energy of the 5j8j5 ELD with the V 6 (5-5-10-5-5) structure shown in Fig. 1j -both having 6 vacancies in the (6 3 6) supercell -and found that the stability of the system can be dramatically enhanced when the vacancies are arranged in a line (3.59 vs. 6.63 eV for ELD and V 6 , respectively). In analog to graphene 35 , ELDs can significantly modify the electronic properties of silicene, which may transform semimetallic silicene into metallic.
Adsorption and diffusion of silicon atoms on silicene. Adatoms are also an important defect of silicene. Consistent with previous studies 23, 36 , our calculations found that Si adatoms prefer to adsorb at the top site of silicene, and forms the so-called ''dumbbell configuration'' 23,37 . Figure 5a shows the possible structures for co-adsorption of two adatoms on a silicene sheet, where the first Si is ''fixed'' at site 1, and the second one could stay at sites A to F. We define the adsorption energy between adatoms and silicene, E ad , as E ad 5 (E tot -E silicene -E si )/N, where E tot , E silicene , and E si denote the total energies of the entire system containing N adatoms, the pristine silicene, and an isolated silicon atom, respectively. The analysis of the E ad and E f values shown in Table I indicates that the adsorption of adatoms on silicene is an exothermic process for all cases, and the systems are energetically more favorable than the pristine silicene when the second adatom locates at sites B to C. However, this is not the case when the two migrating adatoms are adjacent, i.e., the second adatom stays at site A. Table I also shows the computed energy barriers for Si adatoms diffusion on the silicene sheet, by using the LST/QST method 38 . Our calculations revealed that the diffusion of an adatom from A to B is significantly easier than to C (0.75 vs. 1.10 eV). Notably, the diffusion barriers for BRD and CRF are almost identical, but are nearly twice as large as that of the ARB pathway.

Electronic and magnetic properties of adatom structures.
Originating from the changes in localized states in sp 2 -bonded materials, adatoms would produce dispersion less metallic band or cause band-gap opening near the Fermi level. Suffering from the self-interaction errors, the standard DFT calculations are known to significantly underestimate the band gap of semiconductors. Moreover, recent work has demonstrated that semilocal generalizedgradient approximation (GGA) calculations would underestimate the relative stability of a spin-polarized configuration, and lead to an incorrect nonmagnetic ground state 39 . Thereby, we have employed the spin-polarized hybrid functional Heyd-Scuseria-Ernzerhof (HSE) for band structure calculations for the energetically favorable defect structures. As illustrated in Fig. 5b, the pristine silicene is a semimetal with the Dirac-like linear dispersion relationship, and its electronic structures can be significantly perturbed due to the adatom-induced symmetry-breaking effects. Indeed, upon the adsorption of a single Si atom at site 1, the Dirac cone vanishes and several bands near the Fermi level become almost flat, implying the strongly localized charge (see Figs. 5a and 5c). The HSE calculations for the single-adatom structure predict a spin-polarized ground state with a spin moment of 2 m B , consistent with the literature 22,37 . Nevertheless, the band gap from the hybrid functional (0.40 eV) is significantly larger than that from GGA (0.08 eV) 37 . Further adsorption of the second Si atom leads to a considerably reduced band gap at site B (0.18 eV from HSE), a similar band gap at site C (0.36 eV), and a vanishing magnetic moment in both systems (see Figs. 5d and 5e). We also notice that all remaining di-adatom structures are metallic, with a spin moment of 2, 4, 0, and 4 m B at sites A, D, E, and F, respectively. We also checked the co-adsorption of three adatoms at sites 1, B, and F, which found to be the ground state with a spin moment of 2 m B . Considering that adatom-induced structural deformation is significant in silicene, and the strain can be used to effectively tune their electronic band structures 40 , it is thus understandable that the electronic and Comparing the thoroughly investigated C adatoms on graphene [41][42] with Si adatoms on silicene studied in the present work, we see their differences and analogies in many aspects. For example, C adatoms prefer to adsorb to the bridge sites of graphene, whilst Si adatoms stay at the top site of silicene, along with a much larger binding energy in the latter case. On the other hand, the smallest energy barrier for C diffuses on graphene (0.74 eV) is close to that for Si diffuses from site A to B on silicene. Moreover, carbon adatoms also cause important change in electronic and magnetic properties, which can vary the spin moment of graphene from 0.12 to 0.27 m B at different coverages and adsorption sites 42 .

Discussion
In summary, we theoretically studied the structure, energetics, and electronic properties of various defects in silicene, in the form of  vacancy clusters, extended line defects, and di-adatoms. We have demonstrated that the vacancy defects in silicene tend to coalesce into a larger vacancy cluster and significantly stabilize the entire system. The formation energy of vacancy clusters can be dramatically decreased when defects are arranged in a line and formed ELDs, and the formation of the d5d7 and t5t7 ELDs in the silicene sheet is significantly easier than in graphene, due to the buckled silicene requires a lower SW transformation barrier than graphene. We also found that the presence of Si adatoms considerably enhances the stability of the pristine silicene, and more importantly, by changing the arrangement of two silicon atoms one can control the magnetic moment and open up the band gap in the silicene system.

Methods
DFT calculations were carried out using the atomic-centered basis set SIESTA code [43][44][45][46][47][48] . The double-j pluspolarization orbitals (DZP) were adopted as atomic orbital basis sets, and the norm-conserving pseudo-potentials were constructed using the Trouiller and Martins scheme 48 . The charge density was projected on a real space grid with an equivalent cutoff 400 Ry to calculate the self-consistent Hamiltonian matrix elements. Geometry optimizations were performed with a residual force threshold of 0.01 eV/Å . The local spin-density approximation (LSDA) 43 was used for vacancy defect calculations, which has been proven to give reasonable reconstructed structures and defect formation energies for 2D materials [35][36] . We used a (6 3 6) supercell, with a vacuum width of 30 Å , to assess the influence of various local defects in silicene. The configurations of silicene with defects were fully relaxed in terms of cell volume and the atomic coordinates with a conjugate gradient algorithm 49 . However, all angles between lattice vectors were constrained during relaxations. We used a Monkhorst-Pack 50 k-point mesh of (8 3 8 3 1) for geometry optimization and (30 3 30 3 1) for electronic properties analysis. To avoid the interactions between imaging vacancy clusters, we used a larger (8 3 8) supercell in the study of coalescence for systems containing 5-and 6-atom vacancies.