Two-dimensional silica : Structural , mechanical properties , and strain-induced band gap tuning

Two-dimensional silica is of rising interests not only for its practical applications as insulating layers in nanoelectronics, but also as a model material to understand crystals and glasses. In this study, we examine structural and electronic properties of hexagonal and haeckelite phases of silica bilayers by performing first-principles calculations. We find that the corner-sharing SiO4 tetrahedrons in these two phases are locally similar. The robustness and resilience of these tetrahedrons under mechanical perturbation allow effective strain engineering of the electronic structures with band gaps covering a very wide range, from of that for insulators, to wide-, and even narrowgap semiconductors. These findings suggest that the flexible 2D silica holds great promises in developing nanoelectronic devices with strain-tunable performance, and lay the ground for the understanding of crystalline and vitreous phases in 2D, where bilayer silica provides an ideal test-bed. VC 2016 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4939279]


I. INTRODUCTION
Two-dimensional (2D) materials that are single-or fewatomic layers thick have driven intensive research interests recently.The unique structural features and opportunities to engineer the atomic structures at nanoscale allow researchers to explore interesting physics in this reduced dimension.Breakthrough applications have also been envisioned along the development of new materials and devices. 1 Compared to other 2D materials such as graphene, hexagonal boron nitride, and transition metal dichalcogenides (TMDs), 2D silica is interesting because the phase transition between its crystalline bulk phase to the disordered one remains as an open issue in the study of glassy solids. 2 The important roles of silica as insulating layers in semiconductor industry also indicate promising applications of 2D silica in layered heterostructures. 3D silica films grown on metal substrate has recently been achieved, where both monolayer and "bilayer," i.e., two hexagonal silicon layers, are explored. 4In contrast to monolayer films that covalently bonded to the metal substrate, the bilayer could be suspended from the substrate or even supported by graphene. 5The bilayer structure exists in both crystalline phase and vitreous modifications with a corner-sharing SiO 4 network as proposed by Zachariasen in 1932. 6More interestingly, with the help of scanning tunneling and transmission electron microscopies, the crystallinevitreous interface of 2D silica can be identified, and the defect formation, structural deformation, phase transition can be excited and imaged in a controlled manner for detailed investigation. 7,8Evidences from these atomically resolved studies show that pentagons and heptagons are the major types of defects in the vitreous phase, similar as its carbon monolayer analogy, graphene.As a result, the hexagonal and 5/7/7/5 haeckelite phases 9,10 of bilayer silica represent two extremes of the bilayer silica structure.The latter structure was first proposed for the carbon monolayers, 9 which consists of extended pentagon and heptagon defects in an ordered arrangement.As a result, the disordered glassy phase of silica bilayer can be considered as a mixed phase between hexagonal and haeckelite.One of the distinct feature of low-dimensional structures is that strain engineering can be feasibly enabled to tune the material properties such as the electronic structures.Considering the insulating nature of crystalline and vitreous silica in the bulk phase, this approach becomes promising in shrinking the gaps into the range of wide-gap or even narrow-gap semiconductors.Understanding the responses of these material forms under mechanical perturbation is thus of interest in the general scope of studying phase transition in 2D materials.

II. METHODS
In this work, we explore the structural and mechanical properties of bilayer silica by performing first-principles calculations using plane-wave basis set based density functional theory (DFT) methods.The Perdew-Burke-Ernzerhof (PBE) parameterization 11 of the generalized gradient approximation (GGA) is used for the exchange-correlation functional.Additional calculations with van der Waals corrections 12,13 are also carried out to validate the current approach, which show only minor effects on the results presented here.For electronic structure calculations, the hybrid functional Heyd-Scuseria-Ernzerhof hybrid functional (HSE06) 13 is also used considering that DFT calculations with GGA usually underestimate the band gap of semiconductors.Projector augmented wave (PAW) potentials 14 are used for ion-electron interactions.The Vienna Ab-initio Simulation Package a) Email: xuzp@tsinghua.edu.cn0021-8979/2016/119(1)/014301/7/$30.00 V C 2016 AIP Publishing LLC 119, 014301-1 (VASP) is employed for all calculations. 15The atomic structure of crystalline bilayer silica in the hexagonal form is constructed following previous experimental characterization (Fig. 1).The 5/7/7/5 haeckelite structures are modified from hexagonal lattice by changing its in-plane topology. 9,10Unit cells containing 12 and 48 atoms are used for hexagonal and haeckelite phases, respectively.A vacuum layer of 20 A ˚is used to isolate the bilayer.For all results presented in this work, an energy cut-off of 520 eV is used for plane-wave basis sets.Monkhorst-Pack (10 Â 10 Â 1 for hexagonal and 8 Â 6 Â 1 for haeckelite) grid points are used for the Brillouin zone sampling.These settings are verified to achieve a total energy convergence of less than 1 meV/atom.For geometry relaxation, the force on atoms and pressure is converged below 0.01 eV/A ˚and 0.01 GPa, respectively.All structural relaxation calculations are performed using the conjugated gradient algorithm.The methodology is validated by reproducing the bulk modulus of a-quartz as 34.4 GPa, which is close to previous experimental measurements and computational calculations. 16The atomic charges are analyzed following Bader's atom-in-molecule approach. 17

A. Atomic structures
The energy calculated for hexagonal and haeckelite crystals differ by 60 meV per SiO 2 unit, and the hexagonal form is more favorable in energy than the metastable haeckelite.This result is consistent with previous reports. 10The lattice constant of the hexagonal cell is 5.31 A ˚, while for the 5/7/7/5 haeckelite cell, the lengths of the two primitive vectors are 10.00, 12.46 A ˚and the angle between them is 51.47 .The internal structural parameters calculated for both phases are summarized in Table I.For both intralayer and interlayer Si-O bonds, the length r Si-O ¼ 1.63 A ˚is all the same for the two phases.In the hexagonal phase, for all oxygen atoms covalently bonded to two silicon atoms, the Si-O-Si angle bends to 140.9 except for the one bridging two tetrahedrons that lies on a straight line.The O-Si-O angle is 109.5 that indicates the nature of tetrahedral symmetry.Both Si-O-Si and O-Si-O angles spread into a wide distribution in the haeckelite phase, where the pentagons and heptagons are significantly distorted compared to its counterpart of carbon. 9,10The formation energies in the hexagonal and haeckelite phases summarized as a function of the areal density in Fig. S1, 49 where we find that the in-plane areal densities for stress-free structures are almost the same, i.e., 12.21 A ˚2 for each pair of SiO 4 tetrahedrons.

B. Mechanical properties
The mechanical responses of hexagonal bilayer silica are measured by simulating uniaxial and biaxial tensile tests with strain e.We carry out biaxial tests (e ¼ e xx ¼ e yy ) by deforming in-plane dimensions of the unit cell, and relaxing the positions of atoms in it.The stress is then calculated as a function of the biaxial tensile strain, as plotted in Fig. 2(a) for the hexagonal crystal.In calculating the stress, one could define a thickness t of the silica bilayer as the height of two SiO 4 tetrahedrons, that is, the interatomic distance between silicon atoms in the two layers.Alternatively, one can use the 2D stress r 2D that equals to rt to avoid the ambiguous definition of t, where r is the stress following its definition in bulk materials.We notice from the simulation snapshots that the hexagonal lattice deforms affinity till the biaxial strain reaches 0.2 at the peak stress.After that the lattice experiences significant structural distortion and reorganization, which breaks down the hexagonal symmetry and results in asymmetric in-plane stress components r xx and r yy .The structural evolution under strain loading is illustrated in Fig. 3 (see details in the supplementary material). 49The thickness of the bilayer structure keeps almost a constant during the biaxial tension, while the Si-O-Si bond angle is flattened due to the tensile strain.Based on the calculated stress-strain relationship, we estimate the 2D areal stiffness of the hexagonal bilayer silica as 117.8 N/m, which is defined as dr 2D /(dA/A) at the limit of zero strain.Here, A measures the in-plane area of silica bilayer.The tensile strength under biaxial loading is 24.7 N/m at the strain-to-failure of 0.2.Furthermore, charge analysis shows that the atomic charge of silicon keeps increasing with the strain and $1.1e electrons is transferred to it at the breaking point of Si-O bonds, which corresponds to the amount of charge transfer in the bilayer silica.Meanwhile, the oxygen atoms that belong to these broken bonds lose atomic charge of $1.0e.
For uniaxial tensile tests, we deform an orthorhombic supercell of 24 atoms anisotropically, with a doubled size compared to the hexagonal unit cell.The supercell is deformed along the armchair (x) or zigzag (y) chain directions in the hexagonal lattice while the dimension of supercell in the other direction is relaxed to zero pressure.The results are summarized in Fig. 4(a).In both cases, the thickness of bilayer, or the height of SiO 4 tetrahedrons, does not show significant change that is similar as the biaxial case as shown in Fig. 3, while the in-plane Poisson's ratios are $0.5 for both armchair and zigzag directions.For uniaxial tension along the armchair direction, the tensile stiffness is calculated to be 136.3N/m.The structure fails at strain of 0.34 with tensile strength of 35.3 N/m, with in-plane Si-O bonds break, which is much higher than the critical strain at peak stress of other 2D materials such as graphene, 0.2 18 (Table II).Along the zigzag direction, the tensile stiffness is 130.5 N/m.It should be noted that although hexagonal lattice is expected to feature an isotropic tensile stiffness in the basal plane, the calculations here are obtained at finite strain, which breaks down the hexagonal symmetry and leads to the slight difference in the tensile stiffness in the armchair and zigzag directions.The tensile strength and strain corresponding to the peak stress are calculated to be 38.3N/m and 0.4.There are tails with finite amplitude of stress carried by the structure after the peak, which may arise due to the limited size of supercells we used in the calculations.
As a two-dimensional material with single-atom thickness, the silica bilayer is a flexible structure in the out-of-plane direction, and thermal fluctuation could result in large-amplitude bending deformation.The three-atomlayer structure of 2D bilayer silica is expected to result in a much higher bending stiffness j compared to other 2D monolayer materials such as graphene and hexagonal boron nitride, where j $ 1 eV.To gain some insights into its outof-plane mechanical properties, we calculate the phonon dispersion of hexagonal silica bilayer.Density functional perturbation theory (DFPT) calculations using a 4 Â 4 supercell are performed using VASP and the PHONOPY package. 19The settings are validated by reproducing the bending stiffness j ¼ 1.6 eV of graphene.Phonon calculations, reflecting soft modes, provide a criterion to quantify the structural stabilities.Our calculations do not find soft modes in the silica bilayers.From the phonon dispersion plotted in Fig. 4(b) (the full phonon dispersion is provided in Fig. S2), 49 we calculate the group velocities of acoustic branches and then quantify mechanical properties of the bilayer.The results show that the group velocities of longitudinal and transverse phonons are v LA , v TA ¼ 8.9, 6.5 km/s, respectively.These values are about half of the values for graphene, i.e., v LA , v TA ¼ 21.0, 12.8 km/s.The flexural phonon ZA mode shows a quadratic dispersion w 2 ¼ (j/q)q 4 , 20 where q is the in-plane mass density of silica bilayer, and the bending modulus j is determined to be 14.4 eV, which is averaged from the value of 12.5 eV along the C-K path and 16.3 eV along the C-M path.This value is one order higher than that of graphene, 1.6 eV, and other 2D materials (Table II), suggesting that the silica bilayer has a more robust planar structure under thermal or mechanical perturbation.
We now turn to discuss the haeckelite phase of bilayer silica.Biaxial loads are applied to the unit cell and the mechanical responses are summarized in Fig. 2(b) as a function of the biaxial strain.We find that the areal elastic stiffness is 84.4 N/m, which is 28.5% lower than that of the hexagonal phase. 10The brittle fracture process splits into two steps at strain of 0.15 and 0.19, respectively.Two in-plane Si-O bonds break subsequently, indicating a biaxial tensile strength of 20.6 N/m that is 16.8% lower than the hexagonal phase and a stress plateau after the first Si-O bond breaks at 11.0 N/m (see details in the supplementary material). 49niaxial tensile tests are also performed along the directions in perpendicular and parallel to lattice vector a 2 (Fig. 1(b)), which predict tensile stiffness of 84.3 and 114.8 N/m.The tensile strength and strain-to-failure are 29.4,27.6 N/m, and 0.3, 0.24, respectively.In contrast to the hexagonal phase, the haeckelite lattice breaks abruptly after the stress peak.

C. Electronic structures
One of the key properties of silica is its wide band gap.The electronic structures of silica bilayer are thus studied here to explore the effect of its 2D nature using the GGA functional first.Band structure and density of state (DOS) analysis (Figs. S3 and S4) 49 suggests that both stress-free hexagonal and haeckelite crystals are insulators with large band gaps of 5.48 and 5.73 eV, which is close to the value of 5.66 eV for a-silica.These results align with the fact that the local bonding structures are similar in both 2D and bulk silica.Interestingly, we find applying in-plane tensile strain can significantly reduce these wide band gaps.The conduction band minimum (CBM) and valence band maximum (VBM) for both phases of silica bilayers are located at the high-symmetry point resulting in direct band gaps.Under uniaxial or biaxial tensile strain, key features of the band structures remain almost the same, but the valence and conduction bands shift towards each other, followed by an insulator-semiconductor-metal transition at the critical values of strain.The charge distribution of CBM and VBM indicates that changes in the band gaps are mainly due to the competition between near-band-edge states (Fig. S5). 49As summarized in Fig. 5, DFT calculation results using GGA show that under biaxial strain larger than $15%, the band gaps of the hexagonal and haeckelite phases become lower than 2 eV.Uniaxial tensile strain beyond $30% along the armchair and zigzag directions, where the structure is still stable, could reduce the gaps of the hexagonal phase into values lower than that of silicon.Under tension along the zigzag direction, the silica bilayer could even become metallic.HSE06 based calculations, on the other hand, predict $1.2-1.8eV higher band gaps than that from GGA.However, the general trend of gap reduction is similar, and the band gaps are reduced down to the semiconductor range under tensile strain.Compared to other strategies such as doping and applying electric field, this characterized range of strain-induced band gap tuning in 2D silica bilayer is the widest one that has been reported yet (Table III).
Considering the high extensibility of 2D silica as demonstrated in our study, this strain-induced wide-to-narrow transition of band gaps hold great promises in developing 2D nanoelectronics.

D. Additional remarks
Finally, from a practical standpoint, we add some remarks here on the experimental realization of the strain engineering.The defective structures in vitreous bilayer silica observed in experiments could be formed during the growth or external excitation.The mechanical response of bilayer silica explored in this work is in general brittle.However, experimental observation suggests that a 5/7/5/7 cluster in defective bilayer silica can transform into a 6/6/6/6 cluster. 8It should be noted that our first-principles structural relaxation calculations are performed without consideration of thermal fluctuation.At high temperature or under high-energy excitations such as electron irradiation, the bond rearrangement and ductile behavior may be activated, as suggested by previous studies on graphene. 18Exploring these behaviors requires intensive finite-temperature first-principles molecular dynamics simulations and awaits further investigation.It should also be noted that our discussions here are focused on isolated bilayers of silica, without consideration of the substrate effects.It was reported that bilayer silica usually grows on surfaces of noble metals. 21Experimental evidence in combination with DFT calculations provide silica bilayers weakly bound to a metal substrate with the adhesion energy of 513 mJ/m 2 , which is slightly lower than the value of 521 mJ/m 2 of interlayer interaction in graphite. 22Our comparative study of silica bilayer supported on the Ru (0001) FIG. 5.The changes in band gaps of 2D silica bilayer under uniaxial and biaxial strain loading for both the hexagonal and haeckelite phases, calculated from DFT using both GGA and HSE functionals.surface shows that the interfacial electronic coupling is in the regime of charge transfer, which is weak enough to perturb the intrinsic band structures of silica bilayer (Fig. S6). 49s a result, in general, our conclusions made here for the isolated silica bilayers are expected to hold for the supported bilayers as well.

IV. CONCLUSION
To summarize, we examine structures and properties of hexagonal and haeckelite phases of 2D bilayer silica crystals by performing first-principles calculations.We find that the metastable haeckelite phase of bilayer silica is slightly less stable compared to the hexagonal phase, with only 60 meV higher energy per SiO 2 unit.The mixture of these two phases provides a model for the glassy silica bilayer characterized in experiments.Under mechanical perturbation, the cornersharing SiO 4 tetrahedrons are quite robust, and the haeckelite phase with pentagons and heptagons has reduced mechanical resistance compared to the hexagonal phase, which is though still comparable with that of a-quartz with Young's modulus of 97.2 GPa and bulk modulus of 36.4GPa. 23The wide band gap of 2D silica is also close to that of bulk silica, which could, however, be significantly reduced into the narrow band gap range under tensile strain.These understandings could help to understand of the crystalline and vitreous phases and their transition in 2D, where bilayer silica provides an ideal test-bed.

FIG. 1 .
FIG. 1. Atomic structures of (a) hexagonal and (b) haeckelite phases of bilayer silica viewed from top, side, and tilt 15 by x and y of the plane, which show the bridged bi-tetrahedral SiO 4 unit structures.The haeckelite phase is composed of pentagons and heptagons that are found in 2D vitreous silica.

Figure S3 .
Figure S3.Band structures of hexagonal and haeckelite silica bilayer under biaxial tensile strain.

Figure S4 .
Figure S4.Densities of states of hexagonal and haeckelite silica bilayer under biaxial tensile strain.

Figure S5 .
Figure S5.Charge distributions of CBM and VBM in hexagonal silica bilayers, calculated using DFT with GGA functional.

Figure S6 .
Figure S6.Decomposed band structures and densities of states for hexagonal silica bilayer supported by Ru (0001) surface, which shows that the electronic coupling between silica and Ru is in the regime of interfacial charge transfer.

TABLE I .
Structural parameters calculated for the hexagonal and haeckelite phases of 2D bilayer silica.
FIG. 2. Mechanical response of (a) hexagonal and (b) haeckelite bilayer silica under biaxial loadings.The average biaxial stress is used to define the tensile strength and stiffness.

TABLE II .
In-plane tensile stiffness K, critical strain at peak stress the in zigzag (armchair) directions e z c (e a [This article is copyrighted as indicated in the article.Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to ] IP: 0%-6% strain, MoS 2 , self-consistent GW0 (scGW0) calculations.