Electric-field induced penetration of edge states at the interface between monolayer and bilayer graphene

# Electric-field induced penetration of edge states at the interface between monolayer and bilayer graphene

Yasumasa Hasegawa and Mahito Kohmoto Department of Material Science, Graduate School of Material Science, University of Hyogo,
3-2-1 Kouto, Kamigori, Hyogo, 678-1297, Japan
Institute for Solid State Physics, University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8581, Japan
July 7, 2011, revised July 3, 2019
###### Abstract

The edge states in the hybrid system of single-layer and double-layer graphene are studied in the tight-binding model theoretically. The edge states in one side of the interface between single-layer and double-layer graphene are shown to penetrate into the single-layer region when the perpendicular electric field is applied, while they are localized in the double-layer region without electric field. The edge states in another side of the interface are localized in the double-layer region independent of the electric field. This field-induced penetration of the edge states can be applied to switching devices. We also find a new type of the edge states at the boundary between single-layer and the double-layer graphene.

###### pacs:
73.22.Pr, 73.20.-r, 73.40.-c, 81.05.ue

## I Introduction

Recently, single-layer and double-layer graphene have been studied both theoretically and experimentallyCastro Neto et al. (2009), because of the interesting properties such as the Dirac pointsNovoselov et al. (2005); Zhang et al. (2005), anomalous Hall effectNovoselov et al. (2005); Zhang et al. (2005); Gusynin and Sharapov (2005); Hasegawa and Kohmoto (2006), and the edge statesFujita et al. (1996); Nakada et al. (1996); Kohmoto and Hasegawa (2007); Castro et al. (2007, 2008a); Ritter and Lyding (2009); Tao et al. (2011). The double-layer graphene has attracted peculiar interestCastro et al. (2010a) due to a band gap controlled by the electric field, which has been predictedMcCann and Fal’ko (2006); McCann (2006) and observedOhta et al. (2006); Castro et al. (2007); Oostinga et al. (2007); Mak et al. (2009); Kuzmenko et al. (2009); Zhang et al. (2009); Li et al. (2009).

The edge states in the single-layer graphene and double-layer graphene have been studied by many authors. In the single-layer graphene the edge states appear at the zigzag edgesFujita et al. (1996); Nakada et al. (1996) and bearded edges. If the system is anisotropic, the edge states also exist at the armchair edgesKohmoto and Hasegawa (2007). The edge states in the double-layer graphene has been studiedCastro et al. (2007, 2008a, 2008b). The edge states at the interface between single-layer and double-layer graphene have also been studied. Transmission across the boundary has been studied theoretically using the effective-mass approximationNilsson et al. (2007); Nakanishi et al. (2010); González et al. (2010), edge states have been studied theoreticallyCastro et al. (2008b); Hu and Ding (2011), and quantum oscillations have been observed in the interfacePuls et al. (2009). Vacancy-induced localized states in the multilayer graphene has been proposedCastro et al. (2010b).

In this paper we study the edge states in the hybrid system of single-layer and double-layer graphene as shown in Fig. 1. We focus on the edge states localized in the boundary between the single-layer and the double-layer regions. We obtain the new edge states localized in one side of the interface between single-layer and double-layer regions with energy . We show an interesting property that the edge states at the boundary between the single-layer and the double-layer regions penetrate into the single-layer region when the electric field is applied perpendicular to the layers.

## Ii model

We assume zigzag edges in both the first and the second layers. Each layer has two sublattices, which we call , , and , as shown in Fig. 1. The left part and the right part of the single-layer regions have and pairs of and sublattices in the x-direction. In the double-layer region there are quartets of sublattices. The left edge in the left single-layer region has only sublattice, which is labeled as , and the right edge in the right single-layer region ( sublattice) is labeled as . There are two boundaries between single-layer and the double-layer regions. These two boundaries are different from each otherCastro et al. (2008b); González et al. (2010); Nakanishi et al. (2010); Hu and Ding (2011). One of the boundaries has and sublattices and the other has and sublattices. We call these boundaries as -boundary and -boundary, respectively. The position of the -boundary is , and the position of the -boundary is . Note that -boundary and the -boundary always appear as a pair if two boundaries are parallel. We assume that the left boundary is the type and the right boundary is the type.

We adopt the tight-binding model, where the hoppings between the nearest sites in the layer (- and -) are taken to be and the interlayer hoppings between the nearest sites (-) are taken to be . We take into account the energy difference between layers ( and ), which is controlled by the electric field perpendicular to the layers.

We apply the same method which we have used in studying the edge states in the single-layer grapheneKohmoto and Hasegawa (2007). Imposing the periodic boundary conditions in the direction, we can take the wave number as the quantum number. For each the Hamiltonian is written as a matrix, where . The eigenstates () in the Schrödinger equation () are vectors with components (the wave functions for the first layer (, , , and , where ) and the wave functions for the second layer (, , , and , where )).

In the single-layer region (, or ), the Schrödinger equation is written as,

 −2tcosky2ΨB1,n−1−tΨB1,n =(E−ϵ1)ΨA1,n, (1) −2tcosky2ΨA1,n+1−tΨA1,n =(E−ϵ1)ΨB1,n. (2)

The equations in the region of double layer () are given by

 −2tcosky2ΨB1,n−1−tΨB1,n =(E−ϵ1)ΨA1,n, (3) −2tcosky2ΨA1,n+1−tΨA1,n −t⊥ΨA2,n =(E−ϵ1)ΨB1,n, (4) −2tcosky2ΨB2,n−1−tΨB2,n −t⊥ΨB1,n =(E−ϵ2)ΨA2,n, (5) −2tcosky2ΨA2,n+1−tΨA2,n =(E−ϵ2)ΨB2,n. (6)

At the left edge of the single-layer region we obtain the equation to be Eq. (2) with and , since there are no sublattice at the left edge in the single-layer region, i.e.,

 −2tcosky2ΨA1,1=(E−ϵ1)ΨB1,0. (7)

Similarly, we obtain the equation at the right edge of the single-layer region to be Eq. (1) with and ,

 −2tcosky2ΨB1,Le−1=(E−ϵ1)ΨA1,Le, (8)

At the -boundary, the equations are obtained by taking and in Eqs. (3), (4), and (6), since there are no sublattices at the -boundary.

 −2tcosky2ΨB1,Lα−1−tΨB1,Lα =(E−ϵ1)ΨA1,Lα, (9) −2tcosky2ΨA1,Lα+1−tΨA1,Lα =(E−ϵ1)ΨB1,Lα, (10) −2tcosky2ΨA2,Lα+1 =(E−ϵ2)ΨB2,Lα. (11)

The equations at the -boundary are obtained by taking and , in Eqs. (3), (4), and (5). Explicitly, the equations at the -boundary are given by

 −2tcosky2ΨB1,Lβ−1−tΨB1,Lβ= (E−ϵ1)ΨA1,Lβ, (12) −2tcosky2ΨA1,Lβ+1−tΨA1,Lβ −t⊥ΨA2,Lβ= (E−ϵ1)ΨB1,Lβ, (13) −2tcosky2ΨB2,Lβ−1−t⊥ΨB1,Lβ= (E−ϵ2)ΨA2,Lβ. (14)

By taking , we obtain that Eqs. (1) and (2) are two independent equations for and , respectively, in the single-layer regions. When , Eqs. (3) - (6) become two sets of coupled equations for (, ) and (, ) in the double-layer region by taking . However, if , these equations cannot be separated into the independent equations for any . This is the origin of the field-induced penetration of the edge states into the first-layer region at the -boundary, as we will show below.

## Iii strictly localized states at ky=π

Since at , Eqs. (1) - (14) are the equations within the same group of , i.e. the states at are strictly localized at the ellipses or circles in Fig. 1, as in the single-layer grapheneKohmoto and Hasegawa (2007).

The energies at in the single-layer regions are obtained from Eqs (1) and (2), as

 Es,±=±t+ϵ1, (15)

with -fold degeneracy. The energies of the strictly localized states in the double-layer region are obtained as the eigenvalues of the matrix

 Md=⎛⎜ ⎜ ⎜⎝ϵ1−t00−tϵ1−t⊥00−t⊥ϵ2−t00−tϵ2⎞⎟ ⎟ ⎟⎠, (16)

and they are obtained to be

 Ed,±,±=ϵ1+ϵ22 ± ⎷(Δϵ2)2+t2+t2⊥2±√(Δϵ)2t2+t2t2⊥+t4⊥4, (17)

where , with -fold degeneracy.

At the left and the right edges we obtain the energy as

 EL=ER=ϵ1. (18)

At the -boundary we obtain the energies of the strictly localized states as the eigenvalues of the matrix

 Mα=⎛⎜⎝ϵ1−t0−tϵ1000ϵ2⎞⎟⎠, (19)

which are obtained as

 Eα,0=ϵ2, (20)

and

 Eα,±=±t+ϵ1. (21)

The energies of the strictly localized states at the -boundary are obtained as the eigenvalues of the matrix,

 Mβ=⎛⎜⎝ϵ1−t0−tϵ1−t⊥0−t⊥ϵ2⎞⎟⎠. (22)

When and , we obtain the energies as

 Eβ0= ϵ2+t2⊥t2+t2⊥Δϵ +t4t2⊥(t2+t2⊥)4(Δϵ)3+O((Δϵ)5) (23) Eβ±= ±√t2+t2⊥+ϵ2+(2t2+t2⊥)2(t2+t2⊥)Δϵ ±t2⊥(4t2+t2⊥)8(t2+t2⊥)5/2(Δϵ)2+O((Δϵ)3). (24)

The eigenstates with the eigenvalues and are obtained as

 (25)

and

 ⎛⎜ ⎜⎝ΨA1,LβΨB1,LβΨA2,Lβ⎞⎟ ⎟⎠ = ΨA1,Lβ⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝1∓√t2+t2⊥t+t2⊥Δϵ2t(t2+t2⊥)+O((Δϵ)2))t⊥t∓t⊥Δϵt√t2+t2⊥+O((Δϵ)2))⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠, (26)

respectively.

Since , , , , and are different from the energies of the macroscopically degenerate states ( and ), the eigenstates with these energies become the well-defined edge states at , as we will show below.

## Iv edge states without perpendicular electric field

### iv.1 edge states with E=0

First, we study the edge states in the case of no external electric field (). We plot the energy as a function of in Fig. LABEL:figsds, where we take , , , , , and . As shown in the previous section, there are four states which have when and , i.e., , , , and . Two of them are edge states localized at each edge in the single-layer regions ( and ) for , same as the single-layer grapheneKohmoto and Hasegawa (2007). The edge state localized at the left edge of the single layer is given by

 ΨB1,n=(−2cosky2)nΨB1,0, (27)

where and other components of are zero. The edge state localized at the right edge of the single layer is given by

 ΨA1,Le−j =(−2cosky2)jΨA1,Le, (28)

where and other components of are zero.

The other edge states with are localized at the and boundaries. As shown in Appendix A, the edge state localized at the -boundary are obtained as

 (ΨB1,Lα+jΨB2,Lα+j)=⎛⎝0(−2cosky2)jΨB2,Lα⎞⎠, (29)

where and other components of are zero. The edge states at the boundary are obtained as

 (ΨA1,Lβ−jΨA2,Lβ−j)=(−2cosky2)jΨA2,Lβ(−t⊥t(1+j)1), (30)

where and other components of are zero.

These results are consistent with the results obtained in the bilayer edgeCastro et al. (2008a) and the graphite stepsCastro et al. (2008b). The edge state at the -boundary has the finite amplitudes of the wave functions at and sites, while that at the -boundary has the finite amplitude only at the sites.

We plot the square of the absolute value of the wave functions in Figs. LABEL:figwavesdse000 (a), (b) and (c), in which we have taken and in order to lift the degeneracy of the edge states. We plot two edge states together in Fig. LABEL:figwavesdse000 (a), which are localized in the left and the right edges of the single-layer regions. There exist two localized states at each boundary between single-layer and double-layer, as shown in Fig. LABEL:figwavesdse000 (b) and (c).

### iv.2 edge states with E≠0

At the -boundary there exist other edge states, which have energy at . In Fig. LABEL:figwavesdse000 (d) we plot the amplitudes of the wave functions of the edge states at and . Note that if is the eigenstate with energy , is also the eigenstate with energy , when . Therefore, , , and for the edge states with are the same as these with .

As seen in Fig. LABEL:figsds (c) and (d), the edge states at exist at . Although the bulk extended states with the same energy also exist at other values of , the density of states has a peak at that energy due to the edge states. Therefore, the edge states can be observed as a peak in the differential conductance () by the spatially resolving scanning tunneling spectroscopy (STS)Ritter and Lyding (2009); Tao et al. (2011) at the -boundary.

These edge states at can be understood by considering the equations in the single-layer regions Eqs. (1) and (2) (see Appendix B).

The existence of these edge states has not been known before, as far as we know. Although the existence of the edge states at has been suggested as a perfectly reflecting states by Nakanishi et al.Nakanishi et al. (2010), the pure edge states are obtained only at in their paper, since they adopted the effective-mass scheme, which can be used only near the Dirac points.

## V edge states in the presence of perpendicular electric field

In this section we study the edge states in the presence of perpendicular electric field. When the electric field is applied perpendicular to the layers, the potential difference between the first and the second layers, , becomes finite. Even in that case we have the edge states at the left and right edges in the single-layer regions, which are given by , and (the edge states at the left edge) or (the edge states at the right edge). The edge states at the -boundary, which are given by and are also not affected by the perpendicular electric field, since the edge states at the -boundary is localized only in the second layer. As seen in Fig. LABEL:figsdse01 (a) and (b), the states with and exist for even when they are in the upper or lower band.

The edge states at the -boundary, however, is changed drastically by the perpendicular electric field and they are quite different from the edge states in the bilayer grapheneCastro et al. (2008a).

We study the edge states at the -boundary with the energy at in the perpendicular electric field in the similar method given in Appendix B.

We study the single-layer region near the -boundary. We assume that the energy difference between the first and second layers, is smaller than the hopping energy in the plane, . Then we obtain

 |e0|=∣∣∣E−ϵ1t∣∣∣≈∣∣ ∣∣t(ϵ2−ϵ1)t2+t2⊥∣∣ ∣∣<1, (31)

when . The eigenvalues () of the matrix (Eq. (48)) in the right single layer-region () are real when

 |a|<|1−|e0||=1−|e0|, (32)

where as discussed in Appendix B. We expand and (defined in Eqs. (50) and (51)) in , and we obtain

 λ+≈1−e20a>1, (33)
 λ−≈a1−e20<1, (34)

and

 V≈(2−2e0−2e02). (35)

Therefore, the eigenvector of with the eigenvalue is the edge state localized at the -boundary given by

 (ΨA1,Lβ+1+jΨB1,Lβ+1+j) =λj−ΨB1,Lβ+1⎛⎝−2e0√D−a2+e20+11⎞⎠ ≈λj−ΨB1,Lβ+1([]c−e01), (36)

where .

If , we obtain from Eq. (25) and we obtain from Eq. (67). In this case, Eq. (36) shows that the edge states with and is localized only in the double-layer region at the -boundary, if , as shown Fig. LABEL:figwavesdse000(c).

If due to the perpendicular electric field, and become finite for the edge states at the -boundary, resulting in the penetration of the edge states into the single-layer region.

In this way the strictly localized state at the -boundary at () and becomes the localized states which have the finite amplitudes both in the single-layer and double-layer regions when . In Fig. LABEL:figsdse01 we plot the energy as a function of in the case of the finite energy difference between the first layer and the second layer. The energy of the edge states at the -boundary depends on as shown in Fig. LABEL:figsdse01. In Fig. LABEL:figwavesde0013 (a) and (b) we plot the wave functions of the edge states at and for and . In contract to the case of (Fig. LABEL:figwavesdse000 (c)), all components of the wave functions are finite in the double-layer region at the -boundary (). The wave function has finite amplitudes also in the single-layer region (), which means the penetration of the edge state at the -boundary is induced by the electric field. When approaches to , the energy of the edge states (the lower blue line) deviates from , as seen in Fig. LABEL:figsdse01. The localization length of the edge states becomes large as approaches to (see Fig. LABEL:figwavesde0013 (c) and (d)). The effects of the -boundary () and the right edge () on the edge state at -boundary are seen in Fig. LABEL:figwavesde0013 (c) and (d).

We mention the bulk energy gap caused by the perpendicular electric field. Though the bulk energy gap is originated in the double-layer region, the system with finite width of the single-layer regions in both side as shown in Fig.1 has the bulk energy gap, since the the states except for the edge states are extended in both single-layer and double-layer regions. When the widths of the single-layer regions become larger, the bulk energy gap becomes smaller, as shown in Fig. LABEL:figfig6 (a) and (b), which show the -dependences of the energies for the systems with different width. Therefore, the critical value of , at which the edge state at the -boundary has the energy between the bulk energy gap, depends on the width of the system. The edge states at and (Fig. LABEL:figwavesde0013 (c) and (d)) have energies between bulk gap () for the system with and , while the edge states at (Fig. LABEL:figwavesde0013 (a)) has the energy outside of the bulk gap. Since STS probes the local density of statesRitter and Lyding (2009); Tao et al. (2011), the edge states can be observed even if they are outside of the bulk gap. Therefore, the field-induced penetration of the edge state into the single-layer region can be observed by STS.

As shown in Fig. LABEL:figsdse01 (a), Fig. LABEL:figfig6 (a) and (b), the edge states at the -boundary have energies in the bulk energy gap at if . In this case, the edge states are partially filled in the less than half-filled systems. Then electric current flow along the -boundary ( direction) in both single-layer region and double-layer region (see Fig. LABEL:figwavesde0013 (a), (c), and (d)). The edge states are partially filled even when the perpendicular electric field is not applied. However, the contribution of the edge states at the -boundary to the electric current is only in the double-layer region in the absence of the perpendicular electric field, since the edge states are localized only in the double-layer region of the -boundary (see Fig. LABEL:figwavesdse000 (c)). Thus, the electric conductivity along the -boundary in the single-layer region at the -boundary is increased by the perpendicular electric field.

## Vi conclusion

We have studied the edge states in the hybrid system of the single-layer and double-layer graphene. By using the tight-binding model, we obtain the analytic solution of the edge states at the boundary between single-layer and double-layer regions when perpendicular electric field is not applied. We found that the edge states at the -boundary are localized only in the second layer, and the edge states at the -boundary have finite amplitude in both layers in the double-layer region.

We also find the new edge states with , which are localized at the -boundary.

When the perpendicular electric field is applied, the edge states at the -boundary are shown to change drastically. The edge states at the -boundary have finite amplitudes at all sites in both regions of the boundary. The penetration of the edge states induced by the electric field can be observed experimentally by STSRitter and Lyding (2009); Tao et al. (2011). The bulk gap becomes smaller for the systems with the wider width of the single-layer regions. The edge states, however, are possible to be observed by STS even though the energy of the edge states are not located in the middle of the bulk gap because the edge states has larger amplitudes of the wave functions than that for the extended states with the same energy in the region near the boundary.

We propose a simple method to observe the electric-field-induced penetration of the edge states. The conductivity between two terminals placed at the single-layer region at the -boundary will become large when the perpendicular electric field is applied, as the edge states at the -boundary penetrate into the single-layer region. This electric-field-induced penetration of the edge states can be used as electrical devices.

## Appendix A analytical solutions for the edge states with E=0 at the α and β boundaries when ϵ1=ϵ2

The equations in the double-layer region ((3) - (6)) are decoupled into two groups, when and . Eqs. (3) and (5) are written as

 (10t⊥t1)(ΨB1,nΨB2,n)=(−2cosky2)(ΨB1,n−1ΨB2,n−1), (37)

where . From this equation we obtain

 (ΨB1,Lα+jΨB2,Lα+j)=(−2cosky2)j(10−t⊥tj1)(ΨB1,LαΨB2,Lα), (38)

where . Since we obtain

 ΨB1,Lα =(−2cosky2)LαΨB1,0, (39)

from the equations in the single-layer region (Eq. (1)), we can take for and . Therefore, the edge state localized at the -boundary are obtained as

 (ΨB1,Lα+jΨB2,Lα+j)=⎛⎝0(−2cosky2)jΨB2,Lα⎞⎠, (40)

where and other components of are zero.

In the same way we obtain from equations (4) and (6),

 (ΨA1,Lβ−jΨA2,Lβ−j)=(−2cosky2)j(1−t⊥tj01)(ΨA1,LβΨA2,Lβ), (41)

where . As in the edge state in the -boundary, we can take , when and . Then we obtain from Eq. (13)

 ΨA1,Lβ=−t⊥tΨA2,Lβ, (42)

and we obtain the edge state at the -boundary as

 (ΨA1,Lβ−jΨA2,Lβ−j)=(−2cosky2)jΨA2,Lβ(−t⊥t(1+j)1), (43)

where and other components of are zero.

## Appendix B analytical solutions for the edge states with E≈Eβ,± at the β boundary

By replacing by in Eq. (1), we can write the equations in the single-layer regions, Eqs. (1) and (2), as

 (e01a0)(ΨA1,n+1ΨB1,n+1)=(0a1e0)(ΨA1,nΨB1,n), (44)

where

 a =−2cosky2, (45) e0 =E−ϵ1t. (46)

In Eq. (44), we can take or , and . Then we obtain

 (ΨA1,n+1ΨB1,n+1)=T(ΨA1,nΨB1,n), (47)

where is a matrix given by

 T=⎛⎝1ae0a−e0aa2−e20a⎞⎠. (48)

If , the matrix is diagonalized by the matrix as

 V−1TV=(λ+00λ−), (49)

where are the eigenvalues of the matrix ,

 λ±=12a(a2−e20+1±√D), (50)
 V=(√D−a2+e20+1−2e0−2e0√D−a2+e20+1), (51)
 V−1=1C(√D−a2+e20+12e02e0√D−a2+e20+1), (52)
 C=2√D(√D−a2+e20+1), (53)

and

 D=(a+e0+1)(a−e0+1)(a+e0−1)(a−e0−1). (54)

Note that

 λ+λ−=1. (55)

If , and . In this case we obtain the extended states, if the boundary conditions at , , , and are satisfied. On the other hand, if , two eigenvalues of are real and either or is smaller than 1. In this case and can decrease as obtained from Eq. 47, if the boundary conditions at are satisfied by the corresponding eigenstate of . Note that is obtained if and only if

 ||a|−|e0||>1. (56)

Now we examine the edge states at the -boundary with . When and , we obtain

 |a|≪1 (57)

and

 |e0|≈√1+(t⊥t)2>1. (58)

Then the inequality Eq. (56) is satisfied and the edge states can exist. In order to examine the edge states with , we expand , and in when , as

 D≈(e20−1)2−2(e20+1)a2, (59)
 λ+≈−ae20−1, (60)
 λ−≈−e20−1a, (61)

and

 V≈(−2e0)(−e011−e0). (62)

In this case are real and

 |λ+|<1<|λ−|. (63)

The eigenvector of the matrix with the eigenvalue is given by the first column of matrix . Therefore, when and () satisfy the equation,

 ΨA1,Lβ+1ΨB1,Lβ+1=−√D−a2+e20+12e0, (64)

we obtain from Eq. 47

 (ΨA1,Lβ+1+jΨB1,Lβ+1+j) =λj+ΨB1,Lβ+1⎛⎝−√D−a2+e20+12e01⎞⎠ ≈λj+ΨB1,Lβ+1([]c−e01), (65)

where . This state is the edge state localized at the -boundary.

At the -boundary the equations for , , , , and are obtained from Eq. (1) with and Eq. (13),

 (e01a0)(ΨA1,Lβ+1ΨB1,Lβ+1)=(0a01e0t⊥t)⎛⎜ ⎜⎝ΨA1,LβΨB1,LβΨA2,Lβ⎞⎟ ⎟⎠. (66)

This equation is written as

 (67)

This boundary condition as well as the condition that the wave functions should decrease exponentially in the double-layer region in the left part (the double-layer region) of the -boundary can be satisfied by adjusting , , , and . It is indeed possible as seen in Fig. LABEL:figwavesdse000(d). These edge states at are localized at the -boundary and the wave functions have the finite amplitudes in both sides (both single-layer and double-layer regions) of the -boundary as shown in Fig. LABEL:figwavesdse000(d).

The above method can be applied to the edge states at the single-layer region. We consider the edge states at the left boundary in the left single region () as an example. The boundary condition at is given by

 (ΨA1,0ΨB1,0)=(0ΨB1,0). (68)

This state can be the eigenvector of only when and the eigenvalue of this eigenvector is . Therefore, the edge states at the left boundary in the single-layer region exist only when and , as obtained in the previous section. The edge states at the right boundary can be similarly studied by considering the inverse matrix of .

## References

• Castro Neto et al. (2009) A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov,  and A. K. Geim,