Strained graphene/hexagonal boron nitride: A firstprinciples study
Abstract
Due to atomically thin structure, graphene/hexagonal boron nitride (G/hBN) heterostructures are intensively sensitive to the external forces and deformations being applied to their lattice structure. It’s been shown that strain can tune electronic properties of G/hBN. Also, moiré structures driven by misalignment of graphene and hBN layers introduce new features to the electronic behavior of G/hBN. Specifically, moiré patterns can magnify the strain effects. Here, we study the electronic properties of diverse stacking faults of G/hBN with and without the exertion of inplane strain where the strain is imposed on both layers, simultaneously. Utilizing ab initio calculation, we observe strain induced valley drifts and band gap modulation. We find that prior to the direction and the magnitude of the strain, these effects are strongly dependent on the miorientational angle of the layers. Moreover, a substantial increase in the band gap, about seven times larger than the unstrained case, is demonstrated to emerge for small nonequibiaxial strain imposition, typically .
pacs:
I Introduction
Fewlayered heterostructures constructed of 2D materials, consisting of graphene and its isomorphs, such as hexagonal boron nitride, transition metal dichalcogenides, etc., are introduced as an alternative to graphene for band gap emerging and their versatile and prosperous properties Georgiou et al. (2012); Lee et al. (2013); Roy et al. (2013). Since different 2D layer has diverse elastic and electronic properties, the final properties of heterostructures are strongly affected by strain and stacking, and thus can be tuned to fit new functionalities Novoselov and Castro Neto (2012). Graphene/hexagonal boron nitride (G/hBN) has surpassed other heterostructures in research since it can be employed as an approach to control the electronic properties of graphene leaving its favorite features like mobility unchanged Gannett et al. (2011); Xue et al. (2011); Tang et al. (2015). This is due to the fact that hBN as substrate possesses flat surface and imposes less charge inhomogeneity in graphene. G/hBN also has demonstrated signatures of tunability of band gap regarding the application of strain, corrugation and misorientation of layers Giovannetti et al. (2007); Bokdam et al. (2014); SanJose et al. (2014); Cosma et al. (2014); Wallbank et al. (2015). To date, different band gap magnitude for G/hBN have been reported by both theoretical and experimental studies. The band gap predicted by early ab initio study Giovannetti et al. (2007) of band structure of commensurate G/hBN was objected by recent theoretical studies Sachs et al. (2011); Xue et al. (2011). The reason was argued to be the intrinsic strains due to ignorance of the lattice mismatch between graphene and the underlying hBN. In fact, the incommensurability effects are shown to be responsible for the cancellation of quasiparticle band gaps in more realistic systems with the inclusion of lattice mismatch Bokdam et al. (2014); Yankowitz et al. (2016). Yet, interestingly, a nearly band gap at the primary graphene Dirac cone for aligned G/hBN has been observed, confirming the fact that the physics of G/hBN is not completely known Wang et al. (2016). It is believed firstly, that the nature of the method of measurements itself, and secondly, the increment of the mass term as a consequence of substantial height variation, inplane strain, and reduction of interlayer distance, are the main reasons for the observation of a rather large band gap compared to previous studies.
Rotationally misaligned neighboring layers, alongside the lattice mismatch of crystallographic structure between graphene and hBN leading to quasiperiodic structures, moiré pattern Xue et al. (2011); Yankowitz et al. (2012); Woods et al. (2014), has also been studied as a tool for modulation of electronic properties Mishchenko et al. (2014); Wallbank et al. (2015); Artaud et al. (2016). The modulations include renormalization of Fermi velocity, the emergence of Secondary Dirac cones and band gaps. Yankowitz et al. (2012); Hunt et al. (2013); Wang et al. (2016) Also, commensurateincommensurate transition in moiré patterns as a result of lattice adjustment of graphene to its substrate is shown, both theoretically and experimentally Woods et al. (2014); van Wijk et al. (2014); Jung et al. (2015). These transitions can accumulate strains leading to modification of the electronic properties of G/hBNWoods et al. (2014).
In this paper, using ab initio calculation, we study the consequences of inplane strain on commensurate G/hBN with large misorientation angle, when the strain is applied on both layers simultaneously. The main idea is to study the use of moiré pattern in van der Waals heterostructures for magnification of strain effects, such as modification of the gap energy, especially when the twist angle of the layers becomes large and the superlattice periodicity is typically small (e.g. ). We show that the interplay of the few percent magnitude homogeneous strains and moiré pattern in the experimentally applicable system presents large band gap tunability. This feature can be exploited in 2D material based nanoelectronic devices. Devices in which, the maintenance and successful fabrication of heterostructures with a desired rotational angle of layers and ondemand inplane strain via controlling the commensurability degree, is shown to be experimentally practical Kim et al. (2016); Yankowitz et al. (2016).
Ii Methods
In this section, we briefly address the geometric definition of the commensurate moiré structures. A detailed discussion on the derivation of the commensurate moiré superlattice vectors and the angle at which these structures occur can be found in Ref. Lopes dos Santos et al. (2012) and Ref. Mele (2012). We, then, introduce and calculate electronic properties of three different commensurate structures, both before and after applying strain.
For two layers of honeycomb lattices that are rotated with respect to each other around the normal vector to their planes, commensurate structures take place in discrete angles . These angles are in fact, those at which the lattice translation vectors of the upper and lower layers, and , on the span of their primitive vectors, , become equal, i.e. . This yields that starting from fixed sublattices of both layers at origin, the next sublattice of each layer meet when , with denoting the position of the sublattice in cell of upper layer and being the sublattice vector in lower layer. Accordingly, it can be shown that, the total number of disclosed atoms in commensurate supercell is , and the relative misorientation angle is Uchida et al. (2014). Here we study three nonequivalent commensurate structures (See Fig.1), listed below \bibnote[]The 1.8% lattice mismatch between the graphene and hBN layer is been disregarded to ease the simulation of G/hBN according to the current scope of the possibility of ab initio calculations van Wijk et al. (2014):

() commensurate supercell with misalignment angle of containing atoms in total.

() commensurate supercell with twist angle of and total number of atoms in supercell.

() commensurate supercell with misalignment angle of constituted of atoms in total.
For the sake of convenience and simplicity, we name the structures , and , henceforth throughout this study, respectively.
We study the electronic behavior of these G/hBN structures performing firstprinciples calculations implemented in SIESTA code Soler et al. (2002). Double polarized basis (DZP), with Normconserving pseudopotential, alongside the vdW exchangecorrelation functional, within the conjugate gradient method is employed. For each structure demonstrated in Fig.1, lattice constant, kgrid and mesh cutoff convergence tests have been done, to ensure the consumption of best sets of inputs while same optimized basis is used for relaxation procedure. In this regard, MonkhorstPack kpoint grids and mesh cutoff are used for moiré structures , and , respectively. All simulations include a vacuum space of approximately 20 Å height to exclude any interactions between spurious images of the G/hBN. For an unstrained system, both atomic coordinates and lattice vectors are allowed to relax so that the forces on each atom become less than 0.04 eV/Å.
Our DFT computation shows that the mean value of the relaxed interlayer distance between graphene and hBN for , and is 3.5096 Å, 3.4425 Å, and 3.4573 Å, severally. The difference between the interlayer distances can be attributed to the misorientational angle between the adjacent layers. In fact, the diverse stacking configurations result in different long and shortrange interactions between the layers and the vertical relaxation strains of G/hBN Jung et al. (2015). Also, our DFT results are in good agreement with previous studies Giovannetti et al. (2007); Bokdam et al. (2014). Among symmetric stackings of G/hBN (AB, BA, and AA), the AB stacking, i.e. a Nitrogen atom located on top of the center of a graphene hexagon has smaller onsite energy deviations from graphene, and is more energetically favorable van Wijk et al. (2014); Jung et al. (2014). configuration is the most resembling structure to the G/hBN AB stacking since its rotation angle deviates only 2 from that of the G/hBN AB stacking and possesses the largest number of Nitrogen atoms located on graphene hexagon center per supercell. Therefore, it has the smallest interlayer distance among the three.
Strain can affect the electronic properties of G/hBN through the distortion of the lattice structure which leads to the modification of the overlaps of the atomic orbitals. While uniaxial strain distorts the lattice anisotropically, the biaxial strain restores the lattice symmetries and expands the lattice homogeneously along the inplane axes. Moreover, when applying compressive strain, atoms are pushed closer to each other and in extreme cases where the strain is large, corrugations may occur. Here, we study the impacts of inplane strains in both biaxially strained systems and mixed situations where the G/hBN structure is uniaxially strained in one direction and compressed along the other direction. To strain moiré structures within DFT method, we model the lattice vectors for each of the configurations depicted in Fig.1 as , in which is the strain component along the inplane direction, i.e. . Then, we optimize the lattice structure while the lattice vectors are set to be fixed at their strained values and only the atomic coordinates are allowed to relax in accordance with the vectors.
Iii Results and discussions
In low energy regime, where the most striking electronic properties of the graphene emerge, the electronic behavior of the G/hBN is mostly driven from that of graphene since hBN is a large band gap insulator. Furthermore, being an atomically thin semiconductor, graphene is sensitive to the external forces and can be strongly affected by inplane strain being exerted on its lattice structure. Therefore, the electronic properties of the graphene are modified by strain mainly through geometrical changes and modulation of the interatomic distances leading to changes of the overlaps of the atomic orbitals Guinea (2012); Amorim et al. (2016). Also, strain affects the interlayer interactions of G/hBN through expansion and compression of the lattice structure of the individual layers. Here, we calculate strain energy and discuss the G/hBN electronic band dispersion near Fermi energy before and after application of strain. Next, we address gap modulation of G/hBN moiré patterns in accordance to biaxial and nonequibiaxial strains.
Strain energy is defined as the change in total energy of the system after being deformed through the application of strain, i.e. . Persupercell strain energy is, therefore, the strain energy divided by the number of atoms in G/hBN supercell. Persupercell strain energy of G/hBN as a function of applied strain for three commensurate configurations mentioned earlier is illustrated in Fig.2. All three configurations yield same results ranging from zero to for biaxial strain (See Fig.2(a)). Also, the strain energy demonstrates no dependence on the misoriatational angle, as all three configurations show the same trend and similar strain energy for all biaxial strain values. As in the case of biaxial strain, the surface plot of strain energy of G/hBN, displays almost equal response to strain imposition in all configurations where rise in strain energy to occurs when strain approaches (cf. Fig.2(bd)).
Moreover, as it can be seen from surface plots of strain energy, nonequibiaxial strains in which the structure is stretched in one direction and compressed in opposite direction are less efficient in the modulation of the total energy of the system compared to biaxial (compressive) strains.
Fig.2(e) illustrates the strain energy for biaxial and compressive strain being applied on G/hBN for a larger range of strain magnitude. Right (left) inset is the lattice structure of the G/hBN supercell when 12% biaxial (6% compressive) strain is applied. Comparing the effects of biaxial and compressive strain, we see that the strain energy is approximately equal for both strains (cf. strain energy for 6% in Fig.2(e)). Furthermore, while G/hBN experiences severe changes in construction when being subjected to 12% biaxial strain (right inset of panel (e)), the strain energy doesn’t show any signs of entering plasticity region which is expected to occur as in the case of other 2D materials Topsakal et al. (2009); Li (2012); Zhang et al. (2013); Elahi et al. (2015); Wei and Peng (2014). This can be ascribed to the fact that the system under study is not constituted purely of either graphene or hBN. The competition between interlayer interaction to maintain G/hBN structure and the interatomic interaction in graphene to overcome the strain effects imposes great impacts on G/hBN structure. Hence, graphene exhibits no structure maintenance for the strains greater than and decomposes. In contrast, hBN is still stretchable and can endure greater strains so that its total energy still increases in accordance with the applied strain Peng et al. (2012). Therefore, we see an increment in strain energy in general, while the system has already been disintegrated. The compressive strain, on the other hand, imposes a totally opposite effect on the G/hBN. The compression also does enlarge the strain energy as one gradually go beyond the relaxed structure of G/hBN. The exception is that hBN, due to the larger lattice constant and hence having muchextended structure when isolated, resists large contractions. Consequently, G/hBN starts to corrugate when being exposed to biaxial strain (cf. the left inset of Fig.2(e)).
In Fig.3, we study the modifications of electronic properties due to strain and plot the band dispersion near Fermi energy for diverse cases in which the G/hBN is unstrained, equally strained along the two axes () and subjected to nonequibiaxial strain. Also note that the depicted bands are for the strain configurations that are marked with the black triangle, square and circle in the Fig.2(bd).
Prior to changes of the band gap magnitude, we first comprehend the band topology and relocation of the band gap. As discussed earlier, deformation, both biaxial and nonequibiaxial, distorts lattice structure and thus alters the Brillouin zone (BZ). Therefore, the supercell Brillouin zone (sBZ) corners of the unstrained system do not coincide with that of the system under tension. Also, analogous to monolayer graphene Pereira et al. (2009), sBZ corners, named here as points, are no longer equivalent for nonequibiaxialy strained G/hBN. Thus, essentially, because of time reversal symmetry, there are two groups of equivalent sBZ corners. Hence, as illustrated in Fig.3(b), we pick three different paths to cover all high symmetry points in sBZ for band structure calculation. A good example of this feature is depicted in Fig.3(c), where we present the bands for strained configuration (, ) along the paths depicted in Fig.3(b). It can be seen that the minimum difference between valence and conduction bands occurs in the vicinity of the . This relocation of the valley from the sBZ corners is followed by a discontinuity of velocity (kinks) at s. The same behavior is observable in band dispersion of the other configurations, , and when nonequibiaxial strain is exerted on G/hBN (cf. black circle row of the panel (d)). Valley drifts analogous to ones occurring in nonequibiaxial strains are absent in biaxially strained G/hBN as a consequence of the fact that symmetries are not broken (black triangle row of the panel (d)). Nonequibiaxial strains compared to biaxial strains are more efficient in opening band gaps for all G/hBN structures.
For a better understanding of the valley drift, discontinuity of velocity and the relocation of the gap energy, we plot contour maps of the highest valence band, lowest conduction band and the energy difference between them for an exemplary case, in Fig.4. The plots for unstrained moiré structure are brought in the first column and the plots of the second column correspond to strained where and . Strong trigonal symmetry is visible near sBZ corners both for valence and conduction band.
Furthermore, analogous to the monolayer graphene, the valence band maximum and conduction band minimum, move from sBZ corners proportional to the applied strain direction Pereira et al. (2009). As a consequence of displacement of the Dirac cones away from sBZ corners, the gap also drifts away from the and locates beyond the high symmetry path (See Fig.4(f)). Hence, due to the displacement of gaps beyond the high symmetry path, any attempt for observation of strain induced changes of the gap should be made with care.
We find the magnitude of the band gap for the strained G/hBN shown in Fig.4 to be , which in comparison to the unstrained case, (), is times greater. Therefore, gap considerably increases for special cases of inplane strain even though the strain itself is small ( for this case). This is in contrast to the case of external inhomogeneous strains, where it was shown that inhomogeneous strain only has strong effects on DOS and can weakly affect the band gap SanJose et al. (2014). Also, compared to the case of mono and bilayer graphene where spectral gaps emerge for strain magnitudes greater than 20% and only along the preferred direction Pereira et al. (2009); Verberck et al. (2012), moiré structures are more efficient in gap opening and can considerably enhance the strain impacts.
The gap energy as a function of biaxial strain is shown in Fig.5 for , and configurations. The effects of equibiaxial strains on the gap magnitude are typically similar to the case of nonequibiaxial strains except that the anisotropy in the bonds of the uniaxially deformed G/hBN, is absent, leading to relatively smaller changes of the gap energy.
Unlike the resembling behavior of the strain energy for all configurations (Fig.2), band gap energy displays diverse feature for each structure. This implies then, while the strain energy has no rotational dependence, the gap energy is totally affected by the rotation angle between the layers (see Fig.3). Except for the , all configurations represent nonvanishing and relatively small gap energy not exceeding . The gap energy increases rather softer with strain increment, except for configuration. The reason might be due to the variation in the type of sublattice symmetry breaking caused by different misorientation angle alongside the effect of strain on interlayer interactions, similar to the case of asymmetrically strained bilayer graphene Choi et al. (2010). Notice that among the three configurations, the has the largest gap except for the situation where the applied strain is and .
Iv Conclusion
We have shown that exertion of inplane strain on commensurate G/hBN lattice, depending on the direction and the magnitude of the strain and more importantly the misorientational angle results in great impacts on band gap energy. We found that nonsymmetric strains lead to valley drifts and hence relocation of the direct gap from the sBZ corners and the high symmetry path. Therefore, the identification of the true band gap becomes nontrivial and should be done with care. Furthermore, large direct gap energy emerges for small inplane strain imposition, typically for nonbiaxial strains. Provided that the DFT calculation principally underestimates band gaps, band gaps larger than those reported in this study, are expected to be observed in experiments. These achievements offer the use of moiré pattern as a magnifying tool in newly designed nanoelectronic devices where less strain imposition leads to fabrication of these devices experimentally more practical.
V Acknowledgement
Authors acknowledge the support by the Cluster of the school of nano science at Institute for Research in Fundamental Sciences (IPM).
Vi References
References
 Georgiou et al. (2012) T. Georgiou, R. Jalil, B. D. Belle, L. Britnell, R. V. Gorbachev, S. V. Morozov, Y.J. Kim, A. Gholinia, S. J. Haigh, O. Makarovsky, et al., Nat. Nanotechnol. 8, 100 (2012).
 Lee et al. (2013) G.H. Lee, Y.J. Yu, X. Cui, N. Petrone, C.H. Lee, M. S. Choi, D.Y. Lee, C. Lee, W. J. Yoo, K. Watanabe, et al., ACS Nano 7, 7931 (2013).
 Roy et al. (2013) K. Roy, M. Padmanabhan, S. Goswami, T. P. Sai, G. Ramalingam, S. Raghavan, and A. Ghosh, Nat. Nanotechnol. 8, 826 (2013).
 Novoselov and Castro Neto (2012) K. S. Novoselov and a. H. Castro Neto, Phys. Scr. T146, 014006 (2012).
 Gannett et al. (2011) W. Gannett, W. Regan, K. Watanabe, T. Taniguchi, M. F. Crommie, and A. Zettl, Appl. Phys. Lett. 98, 242105 (2011).
 Xue et al. (2011) J. M. Xue, J. SanchezYamagishi, D. Bulmash, P. Jacquod, a. Deshpande, K. Watanabe, T. Taniguchi, P. JarilloHerrero, and B. J. Leroy, Nat. Mater. 10, 282 (2011).
 Tang et al. (2015) S. Tang, H. Wang, H. S. Wang, Q. Sun, X. Zhang, C. Cong, H. Xie, X. Liu, X. Zhou, F. Huang, et al., Nat. Commun. 6, 6499 (2015).
 Giovannetti et al. (2007) G. Giovannetti, P. Khomyakov, G. Brocks, P. Kelly, and J. van den Brink, Phys. Rev. B 76, 073103 (2007).
 Bokdam et al. (2014) M. Bokdam, T. Amlaki, G. Brocks, and P. J. Kelly, Phys. Rev. B 89, 201404 (2014).
 SanJose et al. (2014) P. SanJose, A. GutiérrezRubio, M. Sturla, and F. Guinea, Phys. Rev. B 90, 115152 (2014).
 Cosma et al. (2014) D. a. Cosma, J. R. Wallbank, V. Cheianov, and V. I. Fal’ko, Faraday Discuss. 173, 137 (2014).
 Wallbank et al. (2015) J. R. Wallbank, M. MuchaKruczyÅski, X. Chen, and V. I. Fal’ko, Ann. Phys. 527, 359 (2015).
 Sachs et al. (2011) B. Sachs, T. O. Wehling, M. I. Katsnelson, and A. I. Lichtenstein, Phys. Rev. B 84, 195414 (2011).
 Yankowitz et al. (2016) M. Yankowitz, K. Watanabe, T. Taniguchi, P. SanJose, and B. J. LeRoy, Nat. Commun. 7, 1 (2016).
 Wang et al. (2016) E. Wang, X. Lu, S. Ding, W. Yao, M. Yan, G. Wan, K. Deng, S. Wang, G. Chen, L. Ma, et al., Nat. Phys. 12, 1111 (2016).
 Yankowitz et al. (2012) M. Yankowitz, J. Xue, D. Cormode, J. D. SanchezYamagishi, K. Watanabe, T. Taniguchi, P. JarilloHerrero, P. Jacquod, and B. J. LeRoy, Nat. Phys. 8, 382 (2012).
 Woods et al. (2014) C. R. Woods, L. Britnell, A. Eckmann, R. S. Ma, J. C. Lu, H. M. Guo, X. Lin, G. L. Yu, Y. Cao, R. V. Gorbachev, et al., Nat. Phys. 10, 451 (2014).
 Mishchenko et al. (2014) A. Mishchenko, J. S. Tu, Y. Cao, R. V. Gorbachev, J. R. Wallbank, M. T. Greenaway, V. E. Morozov, S. V. Morozov, M. J. Zhu, S. L. Wong, et al., Nat. Nanotechnol. 9, 808 (2014).
 Artaud et al. (2016) A. Artaud, L. Magaud, T. Le Quang, V. Guisset, P. David, C. Chapelier, and J. Coraux, Sci. Rep. 6, 25670 (2016).
 Hunt et al. (2013) B. Hunt, T. Taniguchi, P. Moon, M. Koshino, and R. C. Ashoori, Science 340, 1427 (2013).
 van Wijk et al. (2014) M. van Wijk, A. Schuring, M. Katsnelson, and A. Fasolino, Phys. Rev. Lett. 113, 135504 (2014).
 Jung et al. (2015) J. Jung, A. M. DaSilva, A. H. MacDonald, and S. Adam, Nat. Commun. 6, 6308 (2015).
 Kim et al. (2016) K. Kim, M. Yankowitz, B. Fallahazad, S. Kang, H. C. P. Movva, S. Huang, S. Larentis, C. M. Corbet, T. Taniguchi, K. Watanabe, et al., Nano Lett. 16, 1989 (2016).
 Lopes dos Santos et al. (2012) J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. B 86, 155449 (2012).
 Mele (2012) E. J. Mele, J. Phys. D: Appl. Phys. 45, 154004 (2012).
 Uchida et al. (2014) K. Uchida, S. Furuya, J.I. Iwata, and A. Oshiyama, Phys. Rev. B 90, 155451 (2014).
 (27) The 1.8% lattice mismatch between the graphene and hBN layer is been disregarded to ease the simulation of G/hBN according to the current scope of the possibility of ab initio calculations van Wijk et al. (2014).
 Soler et al. (2002) J. M. Soler, E. Artacho, J. D. Gale, A. GarcÃa, J. Junquera, P. OrdejÃ³n, and D. SÃ¡nchezPortal, J. Phys. Cond. Matt. 14, 2745 (2002).
 Jung et al. (2014) J. Jung, A. Raoux, Z. Qiao, and A. H. MacDonald, Phys. Rev. B 89, 205414 (2014).
 Guinea (2012) F. Guinea, Solid State Commun. 152, 1437 (2012).
 Amorim et al. (2016) B. Amorim, A. Cortijo, F. de Juan, A. Grushin, F. Guinea, A. GutiÃ©rrezRubio, H. Ochoa, V. Parente, R. RoldÃ¡n, P. SanJose, et al., Phys. Rep. 617, 1 (2016).
 Topsakal et al. (2009) M. Topsakal, E. Aktürk, and S. Ciraci, Phys. Rev. B 79, 115442 (2009).
 Li (2012) T. Li, Phys. Rev. B 85, 235407 (2012).
 Zhang et al. (2013) Q. Zhang, Y. Cheng, L.Y. Gan, and U. Schwingenschlögl, Phys. Rev. B 88, 245447 (2013).
 Elahi et al. (2015) M. Elahi, K. Khaliji, S. M. Tabatabaei, M. Pourfath, and R. Asgari, Phys. Rev. B 91, 115412 (2015).
 Wei and Peng (2014) Q. Wei and X. Peng, Appl. Phys. Lett. 104, 251915 (2014).
 Peng et al. (2012) Q. Peng, W. Ji, and S. De, Comput. Mater. Sci. 56, 11 (2012).
 Pereira et al. (2009) V. M. Pereira, a. H. Castro Neto, and N. M. R. Peres, Phys. Rev. B 80, 045401 (2009).
 Verberck et al. (2012) B. Verberck, B. Partoens, F. M. Peeters, and B. Trauzettel, Phys. Rev. B 85, 125403 (2012).
 Choi et al. (2010) S.M. Choi, S.H. Jhi, and Y.W. Son, Nano Lett. 10, 3486 (2010).