Influence of vacancy defects on the thermal stability of silicene: a reactive molecular dynamics study

Abstract

The effect of vacancy defects on the structural properties and the thermal stability of free standing silicene – a buckled structure of hexagonally arranged silicon atoms – is studied using reactive molecular dynamics simulations. Pristine silicene is found to be stable up to 1500 K, above which the system transits to a three-dimensional amorphous configuration. Vacancy defects result in local structural changes in the system and considerably reduce the thermal stability of silicene: depending on the size of the vacancy defect, the critical temperature decreases by more than 30%. However, the system is still found to be stable well above room temperature within our simulation time of 500 ps. We found that the, stability of silicene can be increased by saturating the dangling bonds at the defect edges by foreign atoms (e.g., hydrogen).

---

1 Introduction

Graphene, a planar sheet of hexagonally arranged carbon atoms, has attracted a lot of attention in recent years due to its unique physical, chemical and mechanical properties.¹ Theory predicts that, similar to graphene, monolayer honeycomb structures can be obtained from the other group-IV elements and their binary compounds.^2,3 Among them the most interesting material is silicene,^4–6 a silicon analogue of graphene, which may be promising for technological applications, in part, due to silicon's compatibility with existing electronics infrastructure. Due to the similarity of the lattice structures, silicene has similar electronic properties as graphene.^6–8 However, unlike graphene, silicene is not stable in a flat configuration: it forms a slightly buckled monolayer structure due to the fact that sp³ hybridization for Si is more stable than sp² hybridization.^4,7,9 On the other hand, such buckling creates new possibilities for controlling the band structure of silicene electrically or by chemical functionalization.^10–16

Since silicene has not been observed in nature and there is no solid phase of silicon similar to graphite, increased efforts have been devoted in recent years to artificially synthesize silicene sheets by implementing sophisticated methods, e.g., chemical exfoliation methods.¹⁷ Another promising approach is to deposit silicon on metal surfaces that interact weakly with the silicon atoms.^18–23 These experiments have provided compelling evidence for the existence of silicene when supported by a substrate, which opens new opportunities for exploring its properties.²⁴ However, for practical applications of silicene, it is important to know its behavior under extreme conditions of temperature²⁵ and strain.^26–29 Here we investigate the effect of vacancy defects on the structural properties and the thermal stability of silicene using reactive molecular dynamics (MD) simulations. Such defects are unavoidable during the fabrication process and strongly affect both the thermal³⁰ and electronic properties of silicene.³¹ On the other hand, defective silicene has recently been proposed as a good candidate for hydrogen separation.³² Such a silicene membrane shows theoretically an enhanced selectivity to hydrogen sorting as compared to its graphene counterpart and is thus promising for gas separation and filtering applications. However, our MD simulations show that such vacancy defects considerably reduce the critical temperature of silicene, above which the system transits to three-dimensional (3D) amorphous structure with dominating sp³ bonding. Depending on the size of the defects, more that 30% reduction of silicene's thermal stability is observed. However, even defective silicene is found to be stable well above room temperature. Moreover, the stability of silicene can be increased if the dangling bonds of the silicon atoms at the defect edges are passivated by, e.g., hydrogen atoms.

2 Theoretical approach

To study the structural (thermal) properties of silicene, MD simulations were performed using the reactive force-field ReaxFF, which, in contrast to classical force-fields, is a general bond-order dependent potential that accounts for bond breaking and bond formation during chemical reactions.³³ The system connectivity is recalculated at every iteration step and non-bounded interactions (van der Waals and Coulomb) are calculated between all atom pairs, irrespective of their connectivity.³³ Since ReaxFF parameters are derived from quantum chemical calculations, it gives energies, transition states, reaction pathways and reactivity trends in agreement with quantum mechanical calculations and experiments.³⁴ Although, the ReaxFF method does not include electrons explicitly, it contains a number of energy contributions (e.g., an over coordination term and a lone-pair term) that directly, albeit implicitly, refer to electrons (see ref. 33 for more details). Furthermore, ReaxFF includes a polarizable charge calculation method, combined with an on-bond Coulomb interaction which enables electrons to move essentially along bonds between atoms with different electronegativity. This reactive force-field has already been successfully used in studying properties of silicon-based nanostructures.³⁵ Our numerical simulations are carried out using the LAMMPS code.^36,37

3 Structural properties

As a typical example, we consider silicene with a supercell containing 960 silicon atoms. We implement periodic boundary conditions along the silicene basal plane to avoid any edge effects. We started our simulations with a flat hexagonal arrangement of silicon atoms – the state which was shown to be unstable in density functional theory (DFT) simulations.⁷ The system transforms to a buckled structure upon geometry optimization as shown in Fig. 1(a) and (b). This buckling is due to reduced π–π* overlap of p_z orbitals, resulting in mixed sp² and sp³ hybridization.⁷ In the optimized geometry the inter-atomic distance is calculated to be d = 2.29 Å, which is close to the recent DFT predictions (d = 2.278 Å).³⁸ The buckling parameter Δ, which shows the out-of-plane displacement of one of the sublattices is Δ = 0.68 Å which is slightly larger than the one obtained in recent DFT reports (Δ = 0.449 Å),³⁸ but is much smaller than the unstable high buckled structure (Δ ≈ 2 Å).⁷

Fig. 1 (a and b) Top and side view of the equilibrium structure of pristine silicene (960 Si atom computational unit cell). (c–f) Optimized structure of silicene with 1 (c), 2 (d), 3 (e) and 4 (f) vacancy defects. Pentagons are highlighted by blue (dark gray) and ring structures with more than 6 silicon atoms are shaded by yellow (light gray).

Next we study structural changes in the system configuration when vacancy defects are present. Such defects can be created by, e.g., using energetic particles such as electrons and/or ions as in the case of other two-dimensional monolayers.³⁹Fig. 2(c)–(f) show the energetically favorable local atomic configurations of the system with up to 4 vacancy defects. It is readily seen from these figures that even a single vacancy defect leads to a significant local distortion of silicon's hexagonal arrangement: a pentagonal ring (highlighted by blue/dark gray) connected with a nine-member ring (yellow/light gray) is formed upon geometry optimization [Fig. 2(c)]. In this 9-member ring only one dangling bond remains, while the unrelaxed structure has three dangling bonds. Another stable configuration for the monovacancy defect is the state with two interconnected pentagonal rings, where all the silicon atoms have three neighbors (not shown here). In the case of a divacancy defect, the system relaxes to the state with an eight-sided elliptical ring which is adjacent to six hexagonal and two pentagonal rings [Fig. 2(d)]. All the atoms in this configuration have threefold coordination. Such a local structural transformation of silicene, which also resembles the local arrangement of carbon atoms in graphene,³⁹ in the presence of monovacancy and divacancy defects has been reported in recent DFT simulations.³¹ Optimized geometries of silicene in the presence of 3 and 4 vacancy defects are shown in Fig. 2(e) and (f). In both cases, two pentagons are formed upon structural optimization, which are connected to a larger ring with 10 and 9 silicon atoms. In the last configuration [Fig. 2(f)] all silicon atoms have threefold coordination [as in the case of the divacancy defect [Fig. 2(d)], while for the 3-vacancy defect [Fig. 2(e)] one silicon atom in the big ring has two neighbors as we also found for the monovacancy defect [Fig. 2(c)]. Thus, the general trend for the reconstructing behavior of silicene is that all the dangling bonds could be healed when an even number of silicon atoms is removed, whereas dangling bonds can remain in the system if the number of vacancies is odd.

Fig. 2 (a) Lindemann index δ of pristine silicene as a function of temperature. Insets show snapshots of the system at T = 1450 K (1) and T = 1550 K (2) after 500 ps simulation time. Silicon atoms with 4 neighbors are shown in red. (b) Time dependence of the averaged bond length d̄ at T = 1550 K. Insets show snapshots of the system at times indicated by the red dots on the d̄(t) curve. (c) Angular distribution of the relative bond orientation at T = 1550 K at the beginning (solid black curve) and at the end of the simulation (dashed red curve).

4 Thermal stability

Starting from the equilibrium state, we increased the temperature of the system up to 2000 K at a rate of 4 K ps⁻¹ using an isothermal-isobaric ensemble with a Nose–Hoover thermostat/barostat for temperature/pressure control. When the desired temperature is reached, constant temperature MD simulations were conducted for 500 ps. The damping constants for temperature and pressure were 100 fs and 2 ps, respectively, and the time step was 0.25 fs in all simulations. For a given temperature we conducted a statistical analysis and the results presented in this manuscript are averaged over an ensemble of 5 different initial distributions of velocities of the atoms. To characterize the thermal stability of the system, we monitored the bond length fluctuations given by the Lindemann index,
where n is the number of Si atoms, δ_i is the Lindemann index of atom i and δ is the Lindemann index for the entire system. During the temperature increase, the system was stable up to T_c = 1650 K, in agreement with recent ab initio molecular dynamics simulations.⁷ However, our 500 ps constant temperature simulations predict slightly smaller critical temperature for silicene, which is clearly seen from Fig. 2(a) where we plotted the Lindemann index δ of the system as a function of temperature. At low temperatures, δ increases linearly with T and the 2D structure of silicene is retained up to T = 1450 K [see inset 1 in Fig. 2(a)]. With further increasing temperature, the system transits to a 3D structure [inset 2 of Fig. 2(a)] with a sharp jump in the δ(T) curve. To characterize this transition, we plotted in Fig. 2(b) the time dependence of the averaged (over all the silicon atoms) bond distance d̄ at T = 1550 K. This process starts with the formation of local “structural defects”, where silicon atoms tend to form local 4-coordinated configurations [red atoms, inset 1 in Fig. 2(b)]. This defective region quickly spreads (within 30 ps) through the rest of the system [inset 2 in Fig. 2(b)] and the system transits to the 3D structure (see inset 2 of Fig. 2(a) and Supplemental online video ref. 40). The average bond distance in the melted structure is d̄ ≈ 2.31 Å, which is close to the experimental bond distance of crystalline silicon d = 2.35 Å.⁴¹ However, the wide angular distribution clearly indicates the amorphous character of the new system [see dashed red curve in Fig. 2(c)]. When we decreased the temperature starting from this 3D structure we did not observe the transition to flat silicene.

In what follows, we study the effect of vacancy defects on the thermal stability of silicene. To our great surprise the stability range of the sample reduced dramatically when introducing a single vacancy defect in our computational unit cell of 960 Si atoms. During our simulation time of 500 ps, the system is found to be stable only up to T_c = 1000 K, which is more than 30% smaller than the critical temperature of pristine silicene [see solid black circles in Fig. 3(a)]. “Structural defects” are formed near the vacancy defect, where the honeycomb arrangement of silicon atoms transforms to a fourfold coordination of silicon atoms [Fig. 3(d)]. These defects expand to the rest of the sample forming defective channels [Fig. 3(e)] that promotes the assembling of silicon atoms into an amorphous 3D structure (similar to the one in the inset 2 of Fig. 2(a); see Supplemental online video ref. 42). Stability range of silicene decreases further with increasing the number of missing silicon atoms (see the results in Fig. 3(a) for n_d = 2–5). Using the criterion of d = 0.05, we calculated the critical temperature T_c at which a structural transformation of the silicene takes place, which is shown in Fig. 3(c). It is seen from this figure that T_c decreases with increasing the size of the defects. In all these cases an additional low temperature plateau is observed in the δ(T) curves due to the formation of local structural defects with fourfold coordination of silicon atoms near the vacancy defects [see Fig. 3(b)]. However, the effect of the vacancies on the thermal stability of silicene can be reduced by saturating the dangling bonds of the silicon atoms near the defects by the hydrogen atoms. For example, for single and divacancy defects we observed an 100 K increase in silicene's stability by hydrogen passivation [compare solid and open circles/squares in Fig. 3(a)].

Fig. 3 (a) Lindemann index δ(T) of silicene for a vacancy defect with different number of missing silicon atoms: n_d = 1 (black circles), n_d = 2 (red square), n_d = 3 (green triangles), n_d = 4 (blue stars) and n_d = 5 (cyan diamonds). Open black circles and open red squares show the results for silicene with hydrogen passivated single and two vacancy defects, respectively. (b) A snapshot of the sample with a divacancy defect at T = 400 K. (c) Melting temperature T_c as a function of n_d constructed using the criterion δ = 0.05. (d and e) Snapshots of the system with monovacancy defects at T = 1000 K and at two different times.

5 Conclusions

We implemented reactive molecular dynamics simulations to study the effect of vacancy defects on the structural properties and the thermal stability of free standing silicene. Vacancy defects result in local structural changes of the silicene lattice, reducing the number of dangling bonds near the defects. They reduce the thermal stability of silicene considerably; more than 30% reduction is observed even for a single vacancy defect in our supercell of 960 Si atoms. When the critical temperature is reached, the quasi-planar silicene transits to a 3D-amorphous configuration with dominating sp³ hybridization. This transition occurs through the formation of “structural defects” in the atomic configuration near the defects with fourfold coordination of the silicon atoms. The stability range of silicene can be increased by passivating the dangling bonds of the silicon atoms near the defects.

Acknowledgements

This work was supported by the Flemish Science Foundation (FWO-Vl) and the Methusalem Foundation of the Flemish Government. The authors are grateful to Prof. Adri van Duin for his support with the ReaxFF force field.

References

1. (a) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva and A. A. Firsov, Science, 2004, 306, 666; (b) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos and A. A. Firsov, Nature, 2005, 438, 197; (c) A. K. Geim and K. S. Novoselov, Nat. Mater., 2007, 6, 183; (d) C. Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov and A. K. Geim, Rev. Mod. Phys., 2009, 81, 109.
2. H. Şahin, S. Cahangirov, M. Topsakal, E. Bekaroglu, E. Akturk, R. T. Senger and S. Ciraci, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 80, 155453.
3. L. Pan, H. J. Liu, Y. W. Wen, X. J. Tan, H. Y. Lv, J. Shi and X. F. Tang, Phys. Lett. A, 2011, 375, 614.
4. K. Takeda and K. Shiraishi, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 14916.
5. Y. Zhang, Y. W. Tan, H. L. Stormer and P. Kim, Nature, 2005, 438, 201.
6. G. Guzmán-Verri and L. C. Lew Yan Voon, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 76, 075131.
7. S. Cahangirov, M. Topsakal, E. Aktuerk, H. Şahin and S. Ciraci, Phys. Rev. Lett., 2009, 102, 236804.
8. S. Lebègue and O. Eriksson, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 115409.
9. U. Röthlisberger, W. Andreoni and M. Parrinello, Phys. Rev. Lett., 1994, 72, 665.
10. Z. Ni, Q. Liu, K. Tang, J. Zheng, J. Zhou, R. Qin, Z. Gao, D. Yu and J. Lu, Nano Lett., 2012, 12, 113.
11. N. D. Drummond, V. Zólyomi and V. I. Fal'ko, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 075423.
12. Y. C. Cheng, Z. Y. Zhu and U. Schwingenschlögl, Europhys. Lett., 2011, 95, 17005.
13. M. J. S. Spencer, T. Morishita, M. Mikami, I. K. Snook, Y. Sugiyama and H. Nakano, Phys. Chem. Chem. Phys., 2011, 13, 15418.
14. N. Gao, W. T. Zheng and Q. Jiang, Phys. Chem. Chem. Phys., 2012, 14, 257.
15. X.-Q. Wang, H.-D. Li and J.-T. Wang, Phys. Chem. Chem. Phys., 2012, 14, 3031.
16. (a) J. Sivek, H. Sahin, B. Partoens and F. M. Peeters, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 085444; (b) H. Sahin and F. M. Peeters, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 085423.
17. Y. Sugiyama, H. Okamoto, T. Mitsuoka, T. Morikawa, K. Nakanishi, T. Ohta and H. Nakano, J. Am. Chem. Soc., 2010, 132, 5946.
18. B. Lalmi, H. Oughaddou, H. Enriquez, A. Kara, S. Vizzini, B. Ealet and B. Aufray, Appl. Phys. Lett., 2010, 97, 223109.
19. C. Leandri, G. Le Lay, B. Aufray, C. Girardeaux, J. Avila, M. E. Davila, M. C. Asensio, C. Ottaviani and A. Cricenti, Surf. Sci., 2005, 574, L9.
20. B. Aufray, A. Vizzini, S. Oughaddou, H. Leandri, C. Ealet and G. Le Lay, Appl. Phys. Lett., 2010, 96, 183102.
21. P. De Padova, C. Quaresima, C. Ottaviani, P. M. Sheverdyaeva, P. Moras, C. Carbone, D. Topwal, B. Olivieri, A. Kara, H. Oughaddou, B. Aufray and G. Le Lay, Appl. Phys. Lett., 2010, 96, 261905.
22. P. Vogt, P. De Padova, C. Quaresima, J. Avila, E. Frantzeskakis, M. C. Asensio, A. Resta, B. Ealet and G. Le Lay, Phys. Rev. Lett., 2012, 108, 155501.
23. H. Jamgotchian, Y. Colignon, N. Hamzaoui, B. Ealet, J. Y. Hoarau, B. Aufray and J. P. Bibérian, J. Phys.: Condens. Matter, 2012, 24, 172001.
24. A. Kara, H. Enriquez, A. P. Seitsonen, L. C. Lew Yan Voon, S. Vizzini, B. Aufray and H. Oughaddou, Surf. Sci. Rep., 2012, 67, 1.
25. T. Morishita, K. Nishio and M. Mikami, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 77, 081401R.
26. M. Topsakal and S. Ciraci, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 024107.
27. G. Liu, M. S. Wu, C. Y. Ouyang and B. Xu, Europhys. Lett., 2012, 99, 17010.
28. A. Iÿnce and S. Erkoc, Phys. Status Solidi B, 2012, 249, 74–81.
29. D. Wu, M. G. Lagally and F. Liu, Phys. Rev. Lett., 2011, 107, 236101.
30. H.-P. Li and R.-Q. Zhang, Europhys. Lett., 2012, 99, 36001.
31. Y.-L. Song, Y. Zhang, J.-M. Zhang, D.-B. Lu and K.-W. Xu, J. Mol. Struct., 2011, 990, 75.
32. W. Hu, X. Wu, Z. Li and J. Yang, Phys. Chem. Chem. Phys., 2013, 15, 5753.
33. A. C. T. van Duin, S. Dasgupta, F. Lorant and W. A. Goddard, J. Phys. Chem. A, 2001, 105, 9396.
34. (a) A. Bagri, C. Mattevi, M. Acik, Y. J. Chabal, M. Chhowalla and V. B. Shenoy, Nat. Chem., 2010, 2, 581; (b) F. Castro-Marcano, A. M. Kamat, M. F. Russo Jr, A. C. T. van Duin and J. P. Mathews, Combust. Flame, 2012, 159, 1272; (c) K. Chenoweth, S. Cheung, A. C. T. van Duin, W. A. Goddard and E. M. Kober, J. Am. Chem. Soc., 2005, 127, 7192; (d) K. Chenoweth, A. C. T. van Duin and W. A. Goddard, J. Phys. Chem. A, 2008, 112, 1040.
35. (a) U. Khalilov, G. Pourtois, A. C. T. van Duin and E. C. Neyts, Chem. Mater., 2012, 24, 2141; (b) U. Khalilov, G. Pourtois, A. C. T. van Duin and E. C. Neyts, J. Phys. Chem., 2012, 116, 8649; (c) U. Khalilov, E. C. Neyts, G. Pourtois and A. C. T. van Duin, J. Phys. Chem., 2011, 115, 24839; (d) E. C. Neyts, U. Khalilov, G. Pourtois and A. C. T. van Duin, J. Phys. Chem., 2011, 115, 4818.
36. G. Kresse and J. Furthmuller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169.
37. S. V. Zybin, W. A. Goddard, P. Xu, A. C. T. van Duin and A. P. Thompson, Appl. Phys. Lett., 2010, 96, 081918.
38. X.-Q. Wang, H.-D. Li and J.-T. Wang, Phys. Chem. Chem. Phys., 2012, 14, 3031.
39. H. Terrones, R. Lv, M. Terrones and M. S. Dresselhaus, Rep. Prog. Phys., 2012, 75, 062501.
40. Supplimental online video: melting process of silicene (960 silicon atoms) at T = 1550 K.
41. N. F. Mott and E. A. Davis, Electronic Processes in Non-Crystalline Materials, Clarendon, Oxford, 1979, p. 323.
42. Supplemental online video: melting process of silicene (960 silicon atoms) with a single vacancy defect in the middle at T = 1000 K.

---