Strain-induced band gaps in bilayer graphene

Strain-induced band gaps in bilayer graphene

B. Verberck, B. Partoens, F.M. Peeters, and B. Trauzettel Departement Fysica, Universiteit Antwerpen, Groenenborgerlaan 171, 2020 Antwerpen, Belgium Institut für Theoretische Physik und Astrophysik, Universität Würzburg, Am Hubland, D-97070 Würzburg, Germany
July 15, 2019

We present a tight-binding investigation of strained bilayer graphene within linear elasticity theory, focusing on the different environments experienced by the A and B carbon atoms of the different sublattices. We find that the inequivalence of the A and B atoms is enhanced by the application of perpendicular strain , which provides a physical mechanism for opening a band gap, most effectively obtained when pulling the two graphene layers apart. In addition, perpendicular strain introduces electron-hole asymmetry and can result in linear electronic dispersion near the K-point. When applying lateral strain to one layer and keeping the other layer fixed, we find the opening of an indirect band gap for small deformations. Our findings suggest experimental means for strain-engineered band gaps in bilayer graphene.

I Introduction

The synthesis of graphene — a single layer of carbon atoms arranged into a honeycomb structure — and measurements of its electronic properties in 2004 by Novoselov et al. Nov1 () has sparked off the development of a whole new research field, continuing to expand rapidly today in both experimental and theoretical directions. The high rate at which the graphene literature is growing is reflected by the regular appearance of reviews on graphene in recent years; for a selection of relevant accounts we refer to Refs. review1 (); Bee (); review2 (); review3 (); review4 (); review5 (); Gei ().

The original excitement about graphene not only came from its two-dimensionality Nov1 (); Nov2 () — flat two-dimensional (2D) crystals had not been successfully fabricated before and were even predicted to be unstable — but also from its electronic properties Nov1 (); Nov2 (); Nov3 (); Kim (). Apart from novel fundamental physics, graphene boasts superior material properties, including high thermal conductivity, current density, carrier mobility, carrier mean free path, strongness, stiffness, elasticity and impermeability.

Not only single graphene sheets but also stacks of a few graphene layers Nov1 () where the layers are coupled by van der Waals interactions as in graphite can be isolated. This has lead to investigations of multilayer graphene, and in particular of bilayer graphene (Fig. 1, top). Interestingly, bilayer graphene displays (almost-)parabolic electronic dispersion at the K-points (Fig. 1, bottom left), making electrons behave differently Kos (); Nil1 (); Par1 (); McC2 () as compared to the single-layer case. Bilayer graphene offers the possibility of applying a bias voltage between the two layers, allowing to tune the band structure. In particular, the inequivalency of the two graphene layers then gives rise to a “Mexican-hat-like” band structure featuring a band gap McC2 (); Gui (); Oht (); McC3 (); Cas () of magnitude (Fig. 1, bottom right), with eV the interplane hopping parameter (see below). A tunable energy gap is important for possible electronic devices.

Figure 1: Graphene bilayer, with A carbon atoms in black and B carbon atoms in gray (top). Electronic spectrum near the K-point in absence (bottom left) and presence (bottom right) of a bias voltage between the layers.

Another means of influencing (mono- or bilayer) graphene’s electronic structure, currently receiving a lot of theoretical attention Gui2 (); Per (); Lee (); Mar (); Muc (); Muc2 (); Nan (); Raz (); Cho (), is provided by mechanical deformations. It has e.g. been shown that it is possible to conceive inhomogeneous strains in single-layer graphene such that they act as a high uniform magnetic field, therefore resulting in strain-induced Landau levels and a zero-field quantum Hall effect Gui2 (). In Ref. Per (), Pereira et al. elaborated a tight-binding description for uniaxially strained single-layer graphene and predicted that a band gap can form upon deformations — along preferred directions — beyond 20%. Recent investigations on strained bilayer graphene include a standard tight-binding treatment of uniaxial strain Lee (), a description of elastic deformations and electron-phonon coupling in bilayer graphene by means of pseudo-magnetic gauge fields Mar (), the effect of strain on the Landau level spectrum and the quantum Hall effect Muc (); Muc2 () and the effect of strain in combination with an external electric field Nan (); Raz ().

In the present paper, we employ a nearest-neighbor tight-binding description and linear elasticity theory to show that a perpendicular strain component modifies the on-site energies of the two carbon sublattices in bilayer graphene, which can open a band gap at the K-point. The band gap is of a different nature than in the case of the “Mexican-hat-like” electronic dispersions. In addition, we consider the case of two graphene layers subjected to different (but uniform) strains, a scenario proposed recently and predicted — by means of ab initio calculations — to also lead to the opening of a band gap Cho (). Uniform deformations parallel to the sheets are studied as well and are compared to the monolayer case Per ().

Ii Strained bilayer-graphene

We first consider the unstrained AB (Bernal) stacking variant of bilayer graphene: two graphene layers (labeled 1 and 2) at Å apart with the A atoms of layer 2 sitting directly on top of the A atoms of layer 1 (Fig. 1, top). The B atoms of layer 1 and the B atoms of layer 2 have no direct neighbor in the opposite layer.

The band structure of bilayer graphene can be described within the tight-binding formalism. The nearest-neighbor tight-binding hamiltonian assumes one free electron provided by each carbon atom and reads


The index stands for the layer (1 or 2), the vectors are the lattice sites of the A atoms of layer . Every A atom in layer 1 is surrounded by three B atoms at relative position vectors , and , where the bond length has a value of Å, and the A atoms in layer 2 have neighboring B atoms at relative positions , . The operators and create and annihilate an electron at site , respectively. The vector ( Å) connects two nearest-neighbor atoms in different graphene sheets. Values for the intra- and inter-plane hopping energies and are obtained from fitting the tight-binding model to experimental data for graphite. Here we use the values quoted in Ref. Chu (): eV and eV. The on-site energies and differ slightly due to the different environments of A and B atoms. The difference eV Chu () is about two orders of magnitude smaller than the hopping parameter and is usually considered only in models going beyond the nearest-neighbor tight-binding hamiltonian (1) where also A–B and B–B hoppings are taken into account (see e.g. Ref. Par1 ()). We therefore put . The value of the on-site energies will turn out to be of critical importance when considering perpendicular strain; we will return to it later. For a comprehensive review on a complete tight-binding description of (unstrained) bilayer graphene, and in particular for a discussion of features in the electronic structure resulting from asymmetry of the diagonal (differences in on-site energies, ), we refer to Ref. Muc3 ().

The direct-space hamiltonian (1) can be converted into a reciprocal-space hamiltonian by introducing four-component spinors


Here, the operators () are the discrete (lattice) Fourier transforms of :


where the -vectors of the first Brillouin zone are defined so that the properties and hold. The hamiltonian matrix reads




The band structure is obtained by solving the secular equation , with the unit matrix, in the first Brillouin zone. At the K and K points, two electron energy bands touch each other at the Fermi level, making the material a semi-metal (zero band gap).

Since the electronic properties of (bilayer) graphene are determined by the band structure near the K point — and —, it is convenient to make a Taylor expansion of around . With and retaining only lowest-order terms in , the hamiltonian matrix near becomes


Here, is defined via . For , the hamiltonian matrix is the complex conjugate of . Note that the quantity can be rewritten as , with the Fermi velocity. Diagonalisation of (6) results in the dispersion shown in the left bottom of Fig. 1.

As mentioned in the introduction, the application of a bias voltage between layers 1 and 2 (replacing the elements and by and and by ) results in a finite band gap McC2 (); Gui (); Oht (); McC3 (); Cas (). However, this band gap lies not at but at a -vector slightly away from (Fig. 1, bottom right), and its theoretical maximum value is . In the following, we will show that applying a strain to the bilayer graphene lattice can also lead to a band gap, which is not bounded by .

In the presence of a displacement field , the position of an atom formerly at is . The effect of a small deformation of the lattice can be written as a correction to the original hamiltonian , with


Within the nearest-neighbor tight-binding approximation, the changes in on-site and hopping energy parameters , , and are related to changes in nearest-neighbor interatomic distances (bond lengths). We stress that it is therefore important to distinguish between A and B sites since they have different environments in the bilayer. As pointed out before, an A site has three in-plane nearest-neighbor B sites and one neighboring A site in the opposite layer at a distance ; a B site has only the three surrounding in-plane A sites as nearest neighbors (see Fig. 1, top). Denoting the three (not necessarily equal) changed bond lengths between neighboring in-plane A and B atoms by (), we have for the corrections to the on-site energies to linear order in the deformations


for each of the two layers. For the hopping parameters we have


In the following we will drop the attributes and . The link between the corrections , , and and the displacement field then comes from considering the bond length corrections and . Details of the calculations and approximations involved are given in Appendix A; the resulting correction to the hamiltonian matrix reads




Here, the quantities () are elements of the strain tensor , the general definition of which involve derivatives of the displacement field components to the coordinates. For the uniform displacements we mostly consider in this work the elements can be related to relative increments/decrements of the bond lengths occuring in the lattice (see Appendix A): .

The leading correction to the hamiltonian matrix at the K-point is -independent and reads


For , the hamiltonian matrix correction is the complex conjugate of . The K-point hamiltonian correction for a strained single graphene layer, of the form


has been derived before in the context of electron-phonon coupling in carbon nanotubes Suz (), where the term proportional to is a deformation potential, a concept going back to Bardeen and Shockley Bar () and the term proportional to corresponds to a bond-length change. However, to our knowledge, the derivation of the strained bilayer hamiltonian [Eqs. (10) – (11)] as given in Appendix A — with particular emphasis on the asymmetry on the diagonal — has not been reported before.

Iii Perpendicular uniform strain

We first consider the bilayer-specific possibility of perpendicular strain, , associated with a change from the inter-layer distance to . Putting the in-plane strain tensor elements to zero, the hamiltonian takes on the form


The phase factors present in Eq. (6) have been eliminated by means of a redefinition of the operators for B sites [Eqs. (58a) – (58b)].

The symmetry-breaking along the diagonal exhibited by the hamiltionian matrix (14) is unusual since it distinguishes not between layers (as e.g. a bias voltage between the two layers would do) but between A and B sublattices. Interestingly, the hamiltonian (14) displays the possibility of the opening of a band gap at : if


the degeneracy of the bands at is lifted, as illustrated in Fig. 2. (Criterion (15), an analytical result, holds when , i.e. when is small compared to .) Recently, Mucha-Kruczyński et al. Muc3 () examined the asymmetry of the (unstrained) graphene bilayer hamiltonian’s diagonal. The possibility of band gap openings and electron-hole asymmetry, as encountered here in Fig. 2, due to different on-site energies, was realized. Here, we show that strain enhances the diagonal’s asymmetry [Eq. (10)] and that elastic deformations therefore provide a physical mechanism for the band structure modifications discussed in Ref. Muc3 ().

Figure 2: Low- band structure for perpendicularly uniformly strained bilayer graphene for various values of .

To check whether the possibility of a strain-induced band gap is experimentally relevant, a more quantitative investigation of criterion (15) is in order. First, the variation of the hopping parameter with interlayer distance can be estimated using Harrison’s relation


which follows from an assumed inverse-square dependence of on Har (). The inequality (15) then becomes


For expansion (), this can be rewritten as


while for contraction (), one obtains


The intra- and inter-plane hopping parameters have values of eV and eV, respectively. The interplane distance is Å, so that eV/Å. Using these values, we have plotted the quantities and in Fig. 3. Interestingly, it follows that the band gap formation criterion is reached more easily (smaller ) in the case of expansion.

Figure 3: (Color online) Plots of the variation of (blue, full line, expansion) and (violet, dashed line, contraction) with . The two horizontal lines represent the values of resulting from using the two extreme values for considered in the present work — eV (olive, dashed-dotted line) and eV (green, dashed line).

To obtain a physically meaningful estimate for the quantity we proceed as follows. First, we recall the tight-binding definition of :


Here, is the electron wave function for a carbon atom, and is the periodic potential of the lattice. Suprisingly, while the tight-binding formalism is a standard method for modelling band structures — especially for carbon structures —, a numerical value for has not yet been established in the literature. The plausible reason is that when the difference between and is neglected (see Sect. II), the presence of on the diagonal of the hamiltonian matrix merely results in an overall shift of the energy bands, making its actual value redundant. In Ref. Rei (), Reich et al. present a detailed comparison of tight-binding and ab initio energy dispersion calculations for graphene and provide fits of the tight-binding parameters to the band structures obtained by density-functional theory. Two variants were considered: one fitting the full -path in -space and one where only -vectors yielding optical transitions with energies less than 4 eV were taken into account, resulting in eV and eV, respectively. We will therefore consider the range eV eV. The simplest way of modelling is to put positively charged ions (charge ) at and :


where spherical coordinates have been introduced. For the carbon electron wave function, we follow the common practice of taking the hydrogen-like orbital:


where Å is the Bohr radius. We then get for the expression


where the proportionality factors left unspecified in Eqs. (21) and (22), together with the factor coming from the azimuthal integration, have been collected into the factor . The integrals


entering Eq. (23) can be solved analytically:


For Å, the integral


has the value Å. From we can then obtain the “calibration value” . For the values eV and eV, we obtain eV/Å and eV/Å, respectively. The dependence of on reads , it is visualised in Fig. 4. Clearly, a linear approximation in the range is justifiable, and the value of can be easily determined numerically. We find eV/Å and eV/Å for eV and eV, respectively; these values are marked by horizontal lines in Fig. 3.

Figure 4: Plots of the variation of with for eV (blue, full line) and eV (violet, dashed line).

Within the present model, it follows that in the case of a large tight-binding on-site energy parameter ( eV), a band gap opens for positive strains larger than (see Fig. 3), corresponding to interlayer distances larger than Å. For , the band gap’s magnitude is about meV (Fig. 2, ).

We recall that formally, criterion (15) is valid when the right-hand side is positive, i.e. when , hence the upper limit of in Fig. 3. Going beyond strains of (both positive and negative) would be well outside the validity of the assumption of small displacements , on which the hamiltonian matrix (10) relies (see Appendix A). For , neglected contributions are of the order which is still acceptable. Similarly, we point out that the unphysical behavior of an ever increasing band gap with increasing must become invalid when linear elasticity, i.e. the assumption of small bond length changes [Eq. (36)], fails.

Based on the foregoing elaborations, stating that the opening of a strain-induced band gap by pulling apart the two graphene layers may be experimentally observed is a fair conclusion. Our main purpose here is not to provide accurate predictions but rather to point out the consequences of the symmetry-breaking along the diagonal in the hamiltonian of strained bilayer graphene. Accurate density-functional theory calculations could provide better numerical estimates, but most importantly, experiments should be undertaken to investigate the effect on the band structure upon pushing together or pulling apart the graphene layers.

Iv Symmetric uniform -strain

We next consider pulling or pushing the bilayer in a direction parallel to the graphene sheets (). The tight-binding hamiltonian then reads


The solutions of the corresponding secular equation are


There can only be a band gap when differs from zero for all -vectors. The same condition arises from the monolayer strain problem, where the solutions of the secular equation read (leading to the Dirac cone at the K-point for ). Pereira et al. Per () have investigated the behavior of in detail. Using the isotropy of a 2D hexagonal lattice, the strain tensor can be written as


with the angle between the tension and the -axis () Perfootnote (), and the Poisson ratio, which takes the graphite value of which we choose in the remainder Bla (). The conclusions made by Pereira et al. are that (i) the minimal strain required for opening a gap is about , that (ii) tension along the zig-zag direction () is optimal for the opening of a band gap, and that (iii) tension along the armchair direction ( never results in a band gap. From Eq. (28) it follows that the inter-plane coupling does not play any role and that the same conclusions are valid for the bilayer.

V Asymmetric uniform -strain

Finally, we consider the possibility of applying different strains (parallel to the bilayer) to the two layers and . In Ref. Cho (), ab initio calculations of a graphene bilayer’s band structure were performed for the situation where one of the two layers is subjected to positive strain along the zig-zag or armchair direction (without the elastic response in the perpendicular direction). Here, we consider a more general situation. We choose to let layer intact (), while layer is pushed or pulled (). Any pushing/pulling direction is allowed; the elastic response is taken into account by using expression (29) for the strain tensor of layer 2.

The deformation of layer introduces a mismatch between the two honeycomb lattices. For infinitesimally small deformations , the crystallographically correct unit cell of the bilayer — the “least common multiple” of the unit cells of the two layers — becomes infinitely large. For practically feasible descriptions of the asymetrically distorted bilayer one is forced to choose deformations that lead to not-too-large unit cells, resulting in a discrete set of strain values . For example, for their ab initio calculations, Choi et al. Cho () considered minimal zig-zag strains of about , requiring a supercell of unit cells in layer and in layer . A tight-binding description faces the same mismatching problem. Apart from the necessity to set up a larger hamiltonian matrix, the inter-layer hopping parameters would have to be carefully reconsidered. Indeed, in the case of asymmetric strain, the A atom originally directly above a particular A atom now has a different position, while a B atom may now have a direct (or almost direct) neighbor above. In a way, the structure becomes a mix of AB and AA stacking Cho (). A similar issue was addressed recently in Ref. San (): spatial modulations in the bilayer interlayer hopping arising due to elastic shearing or twisting were shown to lead to a non-Abelian gauge potential in the description of the low-energy electronic spectrum.

Designing the full tight-binding matrix with correctly interpolated inter-layer hopping parameters is a formidable task, the outcome of which would still only be a limited set of finite accessible strains . As a compromise, we therefore consider the following small- hamiltonian:


where the parameter has to be interpreted as representing the average hopping between the two layers. The application of asymmetric plain affects both the diagonal and the off-diagonal elements; the resulting band structure is therefore non-trivial and worth investigating. Diagonalising the hamiltonian (LABEL:hamIIIa) results in the following secular equation:




Note that the argument does, in general, not cancel out, and that the full 2D dependence of has to be considered. Note also that upon replacing by in Eqs. (31) – (32c), the secular equation remains invariant if changes sign and the phase is shifted by . Hence, the substitution only leads to a band inversion and it suffices, as far as the opening of a band gap concerns, to consider .

The change in band structure comes from the interplay between the values of , and . To obtain an estimate for we use the equivalent of relation (16):


Note that Harrison’s relation [Eqs. (16) and (33)] should be taken as a rule of thumb rather than as an exact result Har (). Other (experimental or theoretical) values than for circulate in the literature (e.g.  Pie (), Bas () or Her ()). As it is our aim to make qualitative conclusions, we have used .

For , we proceed as for in Sect. III and put a charge at and , with any unit vector perpendicular to the -axis (a convenient choice is ) to mimic the periodic lattice potential :


Using the definition (20) for , we numerically calculate the derivative and obtain eV/Å and eV/Å for the cases eV and eV, respectively. Inserting the value eV/Å, we observe the opening of a band gap for arbitrarily small values of . In Fig. 5, low- band structures are shown for a few choices of . The