In search of new reconstructions of (001) α-quartz surface: a first principles study

Centre for Materials Science and Nanotechn Blindern, No.-0316 Oslo, Norway. E-mail: o Singapore University of Technology and De Singapore Department of Physics, University of Oslo, Norway Department of Materials Science and Engin 100 44 Stockholm, Sweden † Electronic supplementary information structures, and calculations for 3 3 surface). See DOI: 10.1039/c4ra10726h Cite this: RSC Adv., 2014, 4, 55599

a-quartz (P3 1 21 space group) is a stable silica (SiO 2 ) polymorph 1 and one of the most abundant minerals, having wide applications in construction, piezoelectric devices, glass production industries, etc. It is also one of the major parts of oil shales 2 and seems to play a signicant role in their properties. A number of studies have been focused on the analysis of a-SiO 2 properties and formed the basic knowledge on phase stability, 1 electronic properties, 3-7 defects, 8-10 some surface properties, 3,[11][12][13][14][15][16] etc. Nevertheless, the understanding of a-SiO 2 surface properties is still weak. It is known that similar to other semiconductors/ insulators, cleaved a-SiO 2 surfaces tend to minimize high energy dangling bonds by surface reconstructions. 16 However, experimental investigations of a-SiO 2 surfaces (especially of reconstructions) are rare and do not provide sufficient details on atomic ordering at the surfaces. It is particularly due to difficulty to use atomic force microscopy (AFM) for the analysis of semiconductor and insulator surfaces. AFM studies reported that a-SiO 2 surface is at. 17 Low-energy electron diffraction (LEED) investigations showed that the a-SiO 2 (001) surface has a (1 Â 1) pattern. 18,19 However, both AFM and LEED studies did not provide any information on the reconstructions at low temperatures. To the best of our knowledge, only ffiffiffiffiffi 84 p Â ffiffiffiffiffi 84 p reconstruction, which takes a place at the temperatures above 873 K, was reported so far. 18 Nevertheless, this reconstruction can be attributed to the a-b quartz phase transition occurring at 846 K. 18 Therefore, the use of computational studies for the understanding of a-SiO 2 surfaces is of signicant interest.
During the last 15 years, a few fundamental studies 11,14,16 showed that a-SiO 2 (001) surface has the lowest energy. In particular, Rignanese et al. reported rst principles study of cleaved and reconstructed a-SiO 2 (001) surfaces, and they suggested that a "dense" surface is the most favorable reconstruction. 14 The properties of the same reconstructed surface were further studied by Rignanese et al. 13 and Goumans et al. 11 Moreover, the "dense" surface was later used for the analysis of molecular adsorption on the SiO 2 surfaces 3,11,13 and SiO 2 interfaces with other materials. 20,21 Recently, using the mix of molecular dynamics (MD) and rst principles calculations, Chen et al. showed that the "dense" surface may further reconstruct. 22 Nevertheless, all previous studies of the reconstructions have some limitations such as small size of the studied systems, limited analysis of system stability, use of empirical interaction potentials for MD, etc. Because of this, in this work, using BOMD simulations and "static" DFT calculations, detailed analysis of the stability of reconstructed and cleaved a-SiO 2 (001) surfaces was performed. We show that the well-known "dense" surface should be considered as metastable. Moreover, the electric and structural properties of the found reconstructed surfaces are analyzed and discussed in details.
All calculations were carried out using the Vienna ab initio simulation package (VASP) with the Perdew-Burke-Ernzerhof (PBE) 23 functional. Projector augmented wave (PAW) pseudopotentials 24,25 were used to model the effect of core electrons. The non-local parts of the pseudopotentials were treated in reciprocal and real space for the "static" DFT and BOMD calculations, respectively. The cutoff energies for the plane wave basis set were set to 400 and 300 eV for all "static" DFT and BOMD calculations, respectively. It should be noted that for BOMD, the increase of the cutoff energy to 400 eV did not affect the nal ("static" DFT) energies of both the 27-layer unit cell and the 2 Â 2 supercell slabs. The Brillouin-zone integrations were performed using G-centered Monkhorst-Pack 26 grids (see Table S1 †).
SiO 2 surface was modelled as a periodic SiO 2 slab containing N (N ¼ 27, 36, or 45) atomic layers and a vacuum region of about 12Å. In this study, we used three slab sizes: unit cell, 2 Â 2 supercell, and 3 Â 3 supercell (only for 27-layer slab, see ESI †). For the studied systems, all atoms were relaxed until the internal forces were smaller than 0.01 eVÅ À1 . The xation of the central layers and its impact on the surface energy of the cleaved surface (for 27-layer slab) were also studied, but the effect was found to be insignicant.
To nd possible reconstructions at a-SiO 2 (001) surface, BOMD simulations with two different annealing/quenching temperature protocols were used. For the rst temperature protocol, atomic velocities were initialized at 1 K using the Maxwell-Boltzmann distribution. Then, the system was heated to 1000 K during the period of 20 ps. Using canonical ensemble BOMD (number of atoms, volume, and temperature are conserved), the resulted system was annealed at 1000 K for a period of 20 ps. Finally, the resulted system was quenched to 1 K with a cooling rate of 0.033 K fs À1 . For the second temperature protocol, the velocities were initialized using the Maxwell-Boltzmann distribution at a few different initial temperatures: 500, 600 (for unit cell calculations only), 1000, and 1500 K. Then, canonical ensemble BOMD simulations were carried out for the period of 15-20 ps with a time step of 1.5-2 fs. To ensure that the increase of simulation time does not affect the structures and energetics, the key BOMD calculations were performed for the period of about 45 ps. Finally, the resulted systems were quenched to 1 K with the cooling rate of 0.033 K fs À1 , and the obtained congurations were further optimized using quasi-Newton algorithm. For all BOMD simulations, the system temperature was controlled using the Nose thermostat. [27][28][29] It should be noted that the results obtained from the rst and the second temperature protocols were found to be similar (the largest energy difference was found to be about 0.02 eV), and hence, in this study, results only for the second temperature protocol are included. Moreover, the key results are presented only for the lowest energy structures predicted based on the mix of BOMD and "static" DFT calculations. Nevertheless, it should be noted that for the slab annealed at different temperatures, the largest difference of "static" DFT energies was found to be as much as 0.02 eV.
In a-SiO 2 , each Si forms bonds with four O atoms with two slightly different bond lengths (in this study, 1.625 and 1.628Å) and the average O-Si-O angle of about 109.5 . The formed Si-O bonds have the covalent polar nature (in some references, 30,31 the Si-O can be also treated as ionic bond). The analysis of computed Bader charges 32 indicates that Si atoms have a charge reduction (electron density reduction) of about 3.23e (the Bader charge is 3.23e), while the O atoms have a charge increase of about 1.62e (the Bader charge is À1.62e). The computed Bader charges correlate well with the analysis of Si and O electronegativities 33 and previously reported results (from 3.26e to 3.33e and from À1.63 to À1.67 for Si and O, respectively). 34,35 For a-SiO 2 , the computed PBE band gap is found to be 5.71 eV, and it is, as expected, smaller compared to experimental values (8.9 AE 0.2 eV for amorphous SiO 2 (ref. 36)), but it is comparable with other theoretical studies (from about 5.6 eV for local density approximation or PBE to 9.4 eV for GW). [5][6][7]14 To study surface stability, the surface energy (g) was calculated as g ¼ (E S À E SiO 2 )/(2A), where E S and E SiO 2 are the energies of the slab and of the bulk SiO 2 , containing the same number of SiO 2 units as the slab; A is the surface area of one side of the slab. Considering the a-SiO 2 structure and taking into account the presence of two surfaces in the slab model, the stoichiometric (001) slab can be built as Si/O-or O/O-terminated slab (see Fig. S1 †). However, since the Si/O-terminated a-SiO 2 (001) slab is highly energetically unstable (the average surface energy was found to be 3.42 J m À2 for the 27-layer slab), this study was limited to the analysis of O/O-terminated slabs. This is consistent with other studies of a-SiO 2 surface stability. 11,14 The surface energy of the cleaved a-SiO 2 (001) surface is found to be 2.23 J m À2 and does not depend on the slab thickness (number of SiO 2 layers). Due to the presence of dangling bonds, the cleaved surface tends to reconstruct. Using BOMD simulation, the evolution of the potential energy of 27layer unit slab (Fig. 1a) was studied. It was found that the potential energy decreases with two steps ((1) reconstruction of one surface and (2) reconstruction of the second surface). As the result of the reconstruction, the under-coordinated Si surface atoms become four-coordinated. Moreover, the comparison of the cleaved and reconstructed structures indicates that the reconstructed system has signicant atomic displacements for the top 6 surface atomic layers (see Fig. 1a). The atomic displacements lead to the formation of 3-and 6-member rings (named based on the number of Si-O pairs in the rings), which are not typical for bulk a-SiO 2 (see Fig. 1a and b). The same reconstruction was also observed for both 36-and 45-layer systems (see ESI †), and for all considered slab thicknesses, the surface energy was found to be 0.39 J m À2 .
The Bader charge analysis shows that atoms at the surface layers of the reconstructed and cleaved surfaces have a signicant difference in the atomic charges. For the reconstructed surface, the average Bader changes for rst-layer O and Si surface atoms are À1.60e and 3.20e, respectively, which are comparable to those (À1.62e and 3.23e) for the bulk a-SiO 2 . In contrast, for the cleaved surface, the average Bader charges on surface O and Si atoms are À1.41e and 3.01e, respectively. This difference correlates well with that for the surface energies of cleaved surface and reconstructed surfaces.
The formed 6-member ring has a triangle-like structure (named as "dense" surface, see Fig. 1b), which is the same as that reported based on MD simulations 11,12 and DFT calculations using the local density approximation. 13,14 Despite the previous studies, the possibility of the "dense" surface reconstruction is not well studied. To the best of our knowledge, only Rignanese et al. 14 reported a rst principles study of surface stability as a function of time, but that study was limited to a short MD simulation period of about 300-400 fs (which is believed to be too short for signicant changes of the reconstruction). Other studies 11,12,22 used a classical MD to obtain the reconstructed surface, which is known to have limitations with describing SiO 2 properties. 37 Because of this, using "dense" prereconstructed 2 Â 2 supercell slab (see Fig. 2a), BOMD simulations of 27-layer slab were performed. For all considered annealing temperatures (500, 1000, and 1500 K), the application of annealing/quenching temperature protocol leads to the transformation of 6-member triangle-like rings to 6-member ellipse-like rings (see Fig. 2a and b, the reconstruction is called as 2 Â 2 reoptimized "dense" surface). This reconstruction is found to be similar to that reported based on classical MD simulations by Chen et al. 22 This reconstruction reduces the surface energy from 0.39 to 0.34 J m À2 for all considered temperatures, indicating that the annealing temperature does not affect the nal energetics (see Fig. 2c). The transformation of 6-member ring induces the turning of the 3-member rings, resulting in atomic displacements at deeper atomic layers than those for the "dense" surface (see Fig. 2a and b). For instance, the analysis of Si-O bond length distributions for the central Si slab layer shows a signicantly broadened distribution compared to that for the "dense" surface (see Fig. 2d). As illustration, for the "dense" 27-layer slab, Si atoms in central  layer form two bonds of mostly identical bond lengths (about 1.628Å; it should be noted that for the "dense" surface, the Si-O bond lengths formed by Si atoms within the central layers approaches to those for the bulk a-SiO 2 (1.625 and 1.628Å) in 45-layer slab). For 2 Â 2 reoptimized "dense" surface in 27-layer model, the Si-O bonds have bond lengths varying from 1.622 to 1.635Å for all considered temperatures (see Fig. 2d). This broadening indicates the existence of indirect surface-surface interaction, which can be induced by overlapping of the atomic displacements. Because of this, it is rough to use the 27layer model for the analysis of reconstruction at real a-SiO 2 (001) surfaces. However, these results can be useful for the analysis of atomic structures of thin a-SiO 2 (001) nanosheets.
Application of annealing/quenching temperature protocol to 36-and 45-layer 2 Â 2 supercell slabs leads to the transformation of 6-member triangle-like rings to 6-member ellipselike rings for all considered temperatures (500, 1000, and 1500 K). The observed rings are structurally similar to those for the 27-layer model. Moreover, similar to the 27-layer model, the formation of the 6-member ellipse-like rings results in the turning of the 3-member rings and signicant atomic displacements within at least 13 (for 45-layer model) atomic layers. Herewith, for the 45-layer model, the largest change of the Si-O bond length for the central layers are about 0.002Å, indicating that the 45-layer slab can be used to represent the real a-SiO 2 (001) surface. For both 36-and 45-layer supercell slabs, the surface energies were found to be 0.35 J m À2 , which are slightly larger compared to that for the 27-layer system (0.34 J m À2 ), but it is still about 0.04 J m À2 lower than that for the "dense" surface (see Fig. 3a). The analysis of Si-O bond length distribution in 6-member rings for the "dense" and 2 Â 2 reoptimized "dense" surfaces suggests that the reconstruction causes the reduction of the dispersion of the distribution (see Fig. 3b and c). As illustration, for the 6-member rings at the "dense" surface, the Si-O bond lengths vary from 1.616 to 1.629 A, while for the 2 Â 2 reoptimized "dense" surface, the variation is from 1.620 to 1.628Å. Moreover, for 2 Â 2 reoptimized "dense" surface, the number of Si-O bonds with bond lengths lying within those for bulk a-SiO 2 (1.625 and 1.628Å) is found to be signicantly larger compared to that for "dense" surface (see Fig. 3b and c). The ring transformation also induces the increase of the average nearest in-plane Si-Si distance (the average nearest in-plane O-O distance is not affected signicantly by the transformation) by about 0.05Å. Taking into account that the Bader charges for the surface atoms at the 2 Â 2 reoptimized "dense" surface are the same as those for the "dense" surface (À1.60e and 3.20e for O and Si, respectively), the increase of Si-Si distance implies the reduction of repulsive electrostatic interaction between the ions. All these are the reasons why the 2 Â 2 reoptimized "dense" surface has a smaller surface energy compared to that for the "dense" surface.
Finally, electronic properties were calculated for the cleaved, "dense", and the 2 Â 2 reoptimized "dense" surfaces (all calculations were done for 45-layer slab supercells, see Fig. 3d). For the cleaved surface, the dangling bonds produce the occupied O-p states above the valence band maximum (VBM) of bulk (d) Density of states for cleaved, "dense", and 2 Â 2 reoptimized "dense" slabs. The Fermi level is set at VBM, and it refers to 0 eV. a-SiO 2 . These states reduce the band gap from 5.71 (for bulk) to 3.68 eV, which is consistent with the observation by Rignanese et al. 14 In contrast, both "dense" and 2 Â 2 reoptimized "dense" surfaces have a similar DOSs and do not provide O-p states above the VBM. The minor difference between the DOSs for the "dense" and 2 Â 2 reoptimized "dense" surface comes from the difference in band gap (5.70 and 5.62 eV for "dense" and 2 Â 2 reoptimized "dense" surface, respectively), indicating that the surface reconstruction can induce the minor reductions of the band gap (by 0.08 eV) due to the surface effects.
In summary, based on DFT calculations, a detailed analysis of possible reconstructions at a-SiO 2 (001) surface was performed. It was found that the "dense" surface has a tendency to reconstruct; the 6-member triangle-like rings at the "dense" surface reconstruct to 6-member ellipse-like rings. This reconstruction is caused by the optimization of both Si-O bond length distribution and Si-Si interactions at the surface layer. The 2 Â 2 reoptimized "dense" surface has the surface energy about 10% lower than the "dense" surface and atomic displacements for the top 13 surface atomic layers. The found lowest energy reconstruction can induce the minor change for the electronic properties (the reduction of band gap by about 0.08 eV).