A Tight-binding model on a honeycomb lattice with a variable lattice spacing

# Zero modes, energy gap, and edge states of anisotropic honeycomb lattice in a magnetic field

## Abstract

We present systematic study of zero modes and gaps by introducing effects of anisotropy of hopping integrals for a tight-binding model on the honeycomb lattice in a magnetic field. The condition for the existence of zero modes is analytically derived. From the condition, it is found that a tiny anisotropy for graphene is sufficient to open a gap around zero energy in a magnetic field. This gap behaves as a non-perturbative and exponential form as a function of the magnetic field. The non-analytic behavior with respect to the magnetic field can be understood as tunneling effects between energy levels around two Dirac zero modes appearing in the honeycomb lattice, and an explicit form of the gap around zero energy is obtained by the WKB method near the merging point of these Dirac zero modes. Effects of the anisotropy for the honeycomb lattices with boundaries are also studied. The condition for the existence of zero energy edge states in a magnetic field is analytically derived. On the basis of the condition, it is recognized that anisotropy of the hopping integrals induces abrupt changes of the number of zero energy edge states, which depend on the shapes of the edges sensitively.

###### pacs:
71.70.Di, 73.43.-f, 81.05.Uw

## I Introduction

Recent experiments on graphene(1); (2); (3); (4); (5) have led to renewed interest in physical properties of electrons on the honeycomb lattice. Despite its simple structure, the honeycomb lattice provides non-trivial physical phenomena which can not be observed in the ordinary square lattice. Among them, much attention has been paid to its peculiar dispersion. In the absence of a magnetic field, the honeycomb lattice has zero modes at the corners K and K’ of the Brillouin zone. By treating these zero modes as Dirac fermions, the unconventional quantization of the Hall conductance observed for graphene was explained(7); (6); (8), although the full proper theoretical treatment of the Hall conductance on the honeycomb lattice was made very recently(9). Moreover, when the system has a boundary, there are edge modes localized on the boundary. The existence of the edge modes depends on a choice of the boundary, and for zigzag and bearded edges there occur large density of states localized on these edges at the Fermi energy due to their flat band structures(10); (11); (12); (13); (14); (15); (16); (17); (18); (19); (20). In addition, it is suggested that the edge modes induce charge accumulation on these edges(21); (22).

In this paper, we study properties of these zero modes in the presence of anisotropy of the hopping integrals in the honeycomb lattice. Recently, the anisotropy of the hopping integrals was introduced by replacing one of the hopping integrals with a general value (23); (24); (25); (26); (9); (28); (27) in order to investigate the unconventional quantum Hall effects on graphene (see Fig.1). (For , we have the anisotropic honeycomb lattice.) In Ref.(9), by using topological arguments, an algebraic expression of the quantum Hall conductance was obtained for almost all gaps including subband gaps, and it was shown that the unconventional quantization of the Hall conductance in a weak magnetic field is realized for weak (), while only the conventional quantization is obtained for strong (). Furthermore, for the graphene case (), the unconventional quantization was found to persist up to the Van Hove singularity(9); (19).

The anisotropy of the hopping parameters is also known to change the peculiar dispersion mentioned above(24); (25); (26); (27). However, in the absence of a magnetic field, its influence is restrictive: Although a gap opens around zero energy for , there remain two zero modes in the Brillouin zone for (24). Therefore, a large anisotropy is needed to change the zero mode structure. As well as the zero mode structure, that of zero energy edge states was shown to change by a large anisotropy(28).

In this paper, it will be shown the situation is drastically changed in the presence of a magnetic field. We analytically derive the condition for the existence of zero modes in a magnetic flux ( and are mutually prime integers), and from the condition it is found that, in the limit of , zero modes exist only for , but a gap around zero energy opens for . In other words, a small anisotropy () is sufficient to open a gap in the presence of a weak magnetic field.

For , the gap around zero energy in a weak magnetic field behaves as a non-perturbative and exponential form as a function of . It will be shown that this behavior is naturally explained in terms of the spontaneous breaking of supersymmetry (29); (30). In particular, an explicit form of the gap around zero energy for is obtained by the WKB method. At , the gap around zero energy in a weak magnetic field is found to make a transition from an exponential (non-perturbative) to a power-law (perturbative) behavior as a function of , and for , the energy bands in a weak magnetic field show linear dependence on .

We will also show that the structure of edge states in the presence of a magnetic field is different from that in the absence of a magnetic field. The condition for the existence of zero energy edge states in a magnetic field is analytically derived, and it is found that the anisotropy of the hopping integrals induces abrupt changes of the number of zero energy edge states, which also sensitively depend on shapes of the edges.

The organization of this paper is as follows. In Sec.II, we present our model. The condition for the existence of zero modes in a magnetic field is analytically derived in Sec.III, both from the secular equation and from the normalizability condition of states with zero energy. On the basis of the condition for the existence of zero modes, the energy spectrum near zero energy in a weak magnetic field is systematically examined in Sec.IV. In Sec.V, zero energy edge states are analyzed, where crucial roles of the anisotropy of the hopping integrals are recognized again. Finally, we summarize our results and discuss possible experimental realization of anisotropy of the hopping integrals in Sec.VI.

## Ii Tight-binding model on the honeycomb lattice in a magnetic field

Let us consider the tight-binding model on the honeycomb lattice with nearest-neighbor hopping in a magnetic field as shown in Fig.1. By denoting wave functions on two sublattices of the honeycomb lattice as and , respectively, the tight-binding model is given by

 Eψn,m = ϕn+1,m−1+e2iπΦnϕn+1,m+1+tϕn,m, Eϕn,m = ψn−1,m+1+e−2iπΦ(n−1)ψn−1,m−1+tψn,m, (1)

where a magnetic flux through the unit hexagon is given by . Here we have introduced anisotropy of the hopping integrals: The hopping integrals of the horizontal bonds are , and those for the other bonds are 1. For simplicity, we neglect the spin degrees of freedom in the following.

## Iii The condition for the existence of zero modes

For the isotropic case (), it was found that zero modes exist for all (rational) values of (31). We now derive the condition for the existence of zero modes in the anisotropic case.

Before examining , let us first consider (24). For , (1) gives

 Eψn,m = ϕn+1,m−1+ϕn+1,m+1+tϕn,m, Eϕn,m = ψn−1,m+1+ψn−1,m−1+tψn,m. (2)

From the Bloch’s theorem, the wave functions are written as 1

 ψn,m=eikxn+ikymψ(k),ϕn,m=eikxn+ikymϕ(k). (3)

Substituting (3) into (2), we have

 Q(k)(ψ(k)ϕ(k))=E(ψ(k)ϕ(k)), (4)

where is given by

 Q(k)=(0D(k)D∗(k)0),D(k)=t+2eikxcosky. (5)

The eigenenergies are given by

 E = ±|D(k)| (6) = ±√(t+2coskxcosky)2+4sin2kxcos2ky.

From (6) with , we find two zero modes at

 k+0:(k0+x,k0+y)=(π,cos−1t2),k−0:(k0−x,k0−y)=(π,−cos−1t2), (7)

for . By expanding and around and in (7),

 kx=k0±x+px,ky=k0±y+py,(|px|,|py|≪1), (8)

is given by

 D±(p)=−itpx±√4−t2py, (9)

where and are those near and , respectively. From (6) and (9), the dispersion relation of the Dirac zero mode is obtained:

 E=±√t2p2x+(4−t2)p2y. (10)

For , the two Dirac zero modes merge into a confluent point

 (k∗x,k∗y)=(π,0), (11)

and for , we have a gap around zero energy.

### iii.1 Derivation from a secular equation

Now we consider . We suppose that is a rational number, ( and are mutually prime integers). Since Eq.(1) has translational symmetry in the -direction, the wave functions are written as

 ψn,m=eikmψn,ϕn,m=eikmϕn, (12)

and (1) becomes

 Eψn = (e−ik+eik+2iπΦn)ϕn+1+tϕn, Eϕn = (eik+e−ik−2iπΦ(n−1))ψn−1+tψn. (13)

By the gauge transformation and , (13) is rewritten as

 Eψn = Anϕn+1+tϕn, Eϕn = A∗n−1ψn−1+tψn, (14)

where with . Since the spectrum is found to be invariant under the transformation , we can restrict the range of to without loss of generality. Moreover, in (14), (,) and (,) obey the same equation, thus from the Bloch’s theorem, we have

 ψn+q=exp(iqθ2)ψn,ϕn+q=exp(iqθ2)ϕn, (15)

where satisfies . Therefore (14) reduces to the eigenequation of a matrix. In the secular equation of this, all non-constant terms containing less than factors of should cancel out each other since the eigenvalue has periodicity with respect to . From this property, it is found that the secular equation is written as the following form:

 F(E2)+f(θ1,θ2)=0, (16)

where is a th-order polynomial of with , and it is independent of . The secular determinant for (14) with determines as

 f(θ1,θ2) = ∣∣ ∣ ∣ ∣ ∣ ∣∣det⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝tA1tA2⋱⋱tAq−1eiqθ2Aqt⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠∣∣ ∣ ∣ ∣ ∣ ∣∣2=∣∣ ∣∣tq+(−1)q−1eiqθ2q∏n=1An∣∣ ∣∣2 (17) = ∣∣∣tq+(−1)q+12cos(q2θ1+q+12π)ei(qθ2+q2θ1+q+12π)∣∣∣2.

Here we have used

 q∏n=1An = q∏n=1[1+ei(θ1+2πΦn)]=1+eiqθ1q∏n=1ei2πΦn (18) = 1+(−1)p(q+1)eiqθ1=1+(−1)q+1eiqθ1,

which is derived from .

When satisfies

 0

the range of is and there exist two independent ’s with . From the secular equation (16), we have two modes at these ’s. On the other hand, if satisfies

 t>21/q, (20)

we have and there is no with . We have a gap around zero energy in this case.

Here we note that the condition for the existence of zero modes for , that is, , is reproduced by (19) with . (When , (1) reduces to that with .)

### iii.2 Derivation from the normalizability condition of states

In Sec. III.1, we derived the condition for the existence of zero modes from the secular equation. Here, we re-derive it from the normalizability condition of states with zero energy.

Let us first consider (14) with ,

 ϕn=−1tAnϕn+1,ψn+1=−1tA∗nψn. (21)

Then for , (18) and (21) lead to

 ϕNq+l=(−1t)q(q∏n=1An)ϕ(N+1)q+l=(−1t)q[1+(−1)q+1eiqθ1]ϕ(N+1)q+l, ψ(N+1)q+l=(−1t)q(q∏n=1A∗n)ψNq+l=(−1t)q[1+(−1)q+1e−iqθ1]ψNq+l, (22)

where . Taking the absolute values of the both sides in (22), we obtain

 |ϕNq+l|=2tq∣∣∣cos(qθ12+q+12π)∣∣∣|ϕ(N+1)q+l|, |ψ(N+1)q+l|=2tq∣∣∣cos(qθ12+q+12π)∣∣∣|ψNq+l|. (23)

For , (23) gives

 |ϕNq+l|<|ϕ(N+1)q+l|,|ψ(N+1)q+l|<|ψNq+l|. (24)

From (24), we see that diverges for , and diverges for . Thus these states are not normalizable, and no relevant zero modes exist. We have a gap around in this case. On the other hand, for , (23) gives

 |ϕNq+l|=|ϕ(N+1)q+l|,|ψ(N+1)q+l|=|ψNq+l|, (25)

at . Thus there exist two zero modes for . These results coincide with those of Sec. III.1.

## Iv Spectrum near zero energy in a weak magnetic field

In this section, we examine the spectrum near zero energy in a weak magnetic field. Although some numerical study was presented in Ref.(23), we perform detailed analytical study here. On the basis of the condition for the existence of zero modes obtained in the previous section, we consider the following four cases separately:

1. , where the condition (19) is always satisfied and we have zero modes for all rational values of .

2. , where zero modes disappear and a gap around opens for .

3. , where one zero mode exists for .

4. , where no zero modes exist.

### iv.1 0<t≤1

We show two examples of the energy bands as a function of in Fig.2. For , we have zero modes for all rational values of . As shown in the following, the energy bands in a weak magnetic field are well described by the continuum approximation.

In the continuum approximation, we use the Landau gauge for a magnetic field . Then, substitution and for (9) with , and (see Appendix A) gives the equation in a weak magnetic field as

 Q±(ψ(x,y)ϕ(x,y))=E(ψ(x,y)ϕ(x,y)), (26)

where is given by

 Q±=(0D±D∗±0),D±=−it^px±√4−t2(^py+πΦx). (27)

Here, and are those near and , respectively. Since and commute each other, we can replace with a c-number . Then putting , we obtain

 Q±(ψ(x)ϕ(x))=E(ψ(x)ϕ(x)), (28)

where

 Q±=(0D±D∗±0),D±=−it^px±√4−t2πΦx. (29)

From (28) and (29), the following equation is obtained:

 H±(ψ(x)ϕ(x))=E2(ψ(x)ϕ(x)), (30)

where

 H±=Q2±=(D±D∗±00D∗±D±). (31)

Therefore, we have

 [t2^p2x+(4−t2)π2Φ2x2−t√4−t2πΦσz](ψ(x)ϕ(x))=E2(ψ(x)ϕ(x)), (32)

around , and

 [t2^p2x+(4−t2)π2Φ2x2+t√4−t2πΦσz](ψ(x)ϕ(x))=E2(ψ(x)ϕ(x)), (33)

around , where is the -component of the Pauli matrix.

Since the equations for in (32) and (33) essentially coincide with those for harmonic oscillators, the energy level for around is given by

 En=±√2πΦt(4−t2)1/4√n,(n=0,1,2,…), (34)

and that around is given by

 En=±√2πΦt(4−t2)1/4√n+1,(n=0,1,2,…). (35)

In a similar manner, the energy levels for around and are given by

 En=±√2πΦt(4−t2)1/4√n+1,(n=0,1,2,…), (36)

and

 En=±√2πΦt(4−t2)1/4√n,(n=0,1,2,…), (37)

respectively. We show the energy levels around and in Fig.4. As illustrated in Fig.3, the energy bands for come to be well fitted by (34) and (35) [or (36) and (37)] in a weak magnetic field.

### iv.2 1<t<2

For , a gap around opens for . This implies that in a weak magnetic field (), a gap around opens by a tiny distortion of graphene, . Since we do not have states, the expressions (34) and (37) need to be modified. We show two examples of energy bands as a function of for in Fig.5.

Let us focus on the gap around . In order to see how it behaves, we plot the natural logarithm of the gap around as a function of in Fig.6. From this, we find that it behaves as

 ΔE∼exp(−α/Φ). (38)

The values of are obtained from Fig.6 as and for and , respectively.

The non-analytic behavior (38) can be understood as breaking of supersymmetry (29); (30); (32) in our model. The operator transforms to and vice versa, which is seen from (28) and (29). By identifying with generators of supersymmetry, in (30) can be considered as sypersymmetric Hamiltonians, ( is “boson”, and is “fermion”). Due to the supersymmetry, there is no perturbative (or power-law) correction with respect to for the states. However, tunneling effects break the supersymmetry spontaneously and the non-perturbative correction (38) appears as a gap around .

When two Dirac zero modes at are close to each other in the momentum space, namely, , the gap around can be estimated by the WKB method. For , the two Dirac zero modes at are located at (7) with

 cos−1t2∼√2−t≡G, (39)

and for , they merge into a confluent point (11). For , it is convenient to expand and around the confluent point (11) instead of or :

 kx=k∗x+px=π+px,ky=k∗y+py=py,(|px|,|py|≪1). (40)

Then in (5) is given by

 D(k)=−2ipx+p2y−G2. (41)

In a weak , we can use the continuum approximation. We use the Landau gauge for a magnetic field . Then, substitution and for (41) with , and (see Appendix A) gives the following equation,

 Q(ψ(x,y)ϕ(x,y))=E(ψ(x,y)ϕ(x,y)), (42)

where is given by

 Q=(0DD∗0),D=−2i^px+(^py+πΦx)2−G2. (43)

In a similar manner as Sec.IV.1, we replace with a c-number and put . Then we obtain

 Q(ψ(x)ϕ(x))=E(ψ(x)ϕ(x)), (44)

where is given by

 Q=(0DD∗0),D=−2i^px+π2Φ2x2−G2. (45)

Identifying with a generator of supersymmetry, we have the supersymmetric Hamiltonian , which satisfies

 H(ψ(x)ϕ(x))=E2(ψ(x)ϕ(x)). (46)

By using the following variable ,

 x+GπΦ=q√1πGΦ, (47)

(46) can be rewritten as

 [−12d2dq2+12q2(1−gq)2−(gq−12)σz](ψ(q)ϕ(q))=E(ψ(q)ϕ(q)), (48)

with

 g=12G√πΦG,E=18πΦGE2. (49)

Therefore, the potential terms for and are given by

 V+(q)=12q2(1−gq)2−gq, (50)

and

 V−(q)=12q2(1−gq)2+gq, (51)

respectively (see Fig.7).

In the leading order of , the potentials (50) and (51) are well approximated by the harmonic oscillator around and . Around , (48) becomes

 [−12d2dq2+12q2+12σz](ψ(q)ϕ(q))=E(ψ(q)ϕ(q)), (52)

and around , (48) becomes

 [−12d2dq2+12(q−1g)2−12σz](ψ(q)ϕ(q))=E(ψ(q)ϕ(q)). (53)

Therefore, the energy levels for around and are given by

 EN+=N+,(N+=0,1,2,…), (54)

and

 EN−=N−+1,(N−=0,1,2,…), (55)

respectively, and those for around and are given by

 EN+=N++1,(N+=0,1,2,…), (56)

and

 EN−=N−,(N−=0,1,2,…), (57)

respectively.

Now take into account tunneling effects between the energy levels around and . The tunneling effects can be estimated by the WKB method presented in Appendix D of Ref.(32). Here we consider only the equation for since the equation for gives the same result. The solution of (52) which vanishes for is given by

where , is a constant, and the parabolic cylinder function(33). The solution of (53) which vanishes for is

where is a constant. We connect these solutions (58) and (59) with that in the forbidden region. In the forbidden region, the usual semi-classical expression for the wave function is available:

 ψ(q)=C1√k(q)exp(−∫qq1k(x)dx)+C2√k(q)exp(∫qq1k(x)dx), (60)

where with in (50), are the turning points, , and are constants. Connecting (58) with (60), and (59) with (60), we obtain

 γ2(−2g2)2E−1Γ(1−E)Γ(−E)=1,γ=e−1/6g2gπ1/2. (61)

For near zero energy (), we have

 (−2g2)2E−1≃−g22,Γ(1−E)≃1,Γ(−E)≃−1E, (62)

thus the solution of (61) for is obtained as

 E(=18πΦGE2)=γ2g22. (63)

This implies that the gap around is given by

 ΔE=4√ΦGexp(−2G33πΦ), (64)

and the exponent in (38) is given by

 α=23πG3=23π(2−t)3/2≡αWKB,(0<2−t≪1). (65)

Now we compare (65) with those obtained numerically for a small . In Fig.8, we show as a function of , and in Table 1, we list them. The relative discrepancy

 δ≡∣∣∣αWKB−αα∣∣∣, (66)

decreases as approaches . This is because the neglected terms in (41) come to be smaller and smaller as approaches .

### iv.3 t=2

At , the two Dirac zero modes at merge into the confluent point (11). As a consequence, a gap around in a weak magnetic field makes a transition from an exponential (non-perturbative) to a power-law (perturbative) behavior as a function of .

In Fig.9(a), we show the energy bands as a function of for . We show the three lowest bands in in the log-log scale for weak magnetic field in Fig.9(b). We fit our data by

 E∼Φκ. (67)

For , we obtain the exponent as , , and for the lowest, the second lowest, and the third lowest bands in , respectively. Thus for we have a power-law behavior .

The behavior is derived analytically from a particular dispersion relation at for (25). For , (6) and (41) with give

 E=±√4p2x+p4y, (68)

which is linear in one direction and quadratic in the other. The exponent is obtained from , where is the area surrounded by an orbit of energy in the momentum space (25). From (68), we have

 S(E)=Γ(1/4)23√2π|E|3/2, (69)

thus .

### iv.4 t>2

For , we do not have zero modes but have a gap around . In Fig.10, we show two examples of the energy bands as a function of for .

Let us study behavior of energy bands in a weak magnetic field by the continuum approximation. For , we expand and around ,

 kx=π+px,ky=py,(|px|,|py|≪1), (70)

then (6) gives

 E=±√E2g+2tp2x+2(t−2)p2y,Eg=t−2. (71)

Thus for small and (), (71) is written as

 E=±(Eg+tt−2p2x+p2y), (72)

which is quadratic in both and . In the continuum approximation, we put and with , and . Then we have

 E=±Eg[1+2π√t(n+12)Φ(t−2)3/2],(n=0,1,2,…), (73)

where we have neglected higher order corrections of . The energy bands in the weak magnetic field limit