Robust Zero Modes in Disordered Two-Dimensional Honeycomb Lattice with Kekulé Bond Ordering
Tohru Kawarabayashi, Yuya Inoue, Ryo Itagaki, Yasuhiro Hatsugai, Hideo Aoki
aa r X i v : . [ c ond - m a t . d i s - nn ] F e b Robust Zero Modes in Disordered Two-Dimensional HoneycombLattice with Kekul´e Bond Ordering
Tohru Kawarabayashi a , Yuya Inoue a , Ryo Itagaki a , Yasuhiro Hatsugai b , Hideo Aoki c,d a Department of Physics, Toho University, Funabashi 274-8510, Japan b Department of Physics, University of Tsukuba, Tsukuba 305-8571, Japan c Department of Physics, University of Tokyo, Hongo, Tokyo 113-0033, Japan d Electronics and Photonics Research Institute, Advanced Industrial Science and Technology (AIST), Tsukuba, Ibaraki305-8568, Japan
Abstract
Robustness of zero-modes of two-dimensional Dirac fermions is examined numerically for thehoneycomb lattice in the presence of Kekul´e bond ordering. The split n = Keywords:
Dirac fermions, chiral symmetry, Kekul´e bond ordering
1. Introduction
Graphene [1, 2] has kicked o ff the physics of Dirac fermions in two dimensions as one ofthe central issues in the condensed-matter physics [3, 4, 5, 6, 7]. In particular, the characteristicquantum Hall e ff ect is a hallmark of the graphene physics. In the developments of theoreticaland experimental understanding of the electronic states for Dirac fermions, there has been a con-siderable progress in the theoretical studies[8, 9, 10, 11] of the charge fractionalization in twodimensions where fractional and irrational charges emerge in association with the topologicaldefects (vortices) of the Kekul´e bond ordering in the honeycomb lattice. Such a system is gen-erally known as a fermion-vortex system having topological zero-energy state in a bulk gap `a laJackiw and Rossi [12]. Experimentally, the Kekul´e bond ordering is now realized not only inmolecular graphene[13] but also in photonic crystals [14] in a controlled manner enabling us toexamine their topological properties.With such a recent progress in mind, here we explore how disorder will a ff ect the topologicalzero-energy state associated with the vortex structure in the Kekul´e bond ordering in two dimen-sions, along with the n = Preprint submitted to Annals of Physics February 2, 2021 . Kekul´e bond ordering
We consider non-interacting electrons on the two-dimensional honeycomb lattice with thenearest-neighbor hopping t . The energy dispersion is given by E ( k ) = ±| d ( ( k ) | with d ( k ) = t [1 + exp( − ik ) + exp( − ik )], where k i ( i = ,
2) denotes the projection k · e i of the wave vectoron the primitive vectors e = ( √ / , / a and e = ( − √ / , / a (Fig. 1(a)) with a beingthe nearest-neighbor distance in the honeycomb lattice. The valence and conduction bands touchwith each other at K and K’ points, ( k , k ) = ( − π/ , π/
3) and (2 π/ , − π/ d ( k ) vanishes and a massless Dirac dispersion emerges. The Kekul´ebond ordering shown in Fig. 1 (b) folds the K and K’ points onto Γ point (0 ,
0) in the Brillouinzone, at which an energy gap ± ∆ K at Γ opens. ∆ K is determined by the degree of the Kekul´ebond ordering (e.g., the di ff erence in the transfer energy between “thick” and “thin” bonds in atight-binding model as depicted in Fig.1). The massless Dirac fermions then acquire a mass (inzero magnetic field), despite the fact that the bond ordering does not break the chiral symmetry.The Kekul´e bond ordering for the honeycomb lattice has been discussed as an instability ina strong magnetic field, where the n = ± ∆ K [17, 18, 19]. Roles of thetopological defects in the Kekul´e bond ordering have also been discussed for transport propertiesof graphene in magnetic fields [20].The tight-binding Hamiltonian for the honeycomb lattice with the nearest-neighbor hoppingis given as H = X < i , j > t i , j e − π i θ i , j c † i c j + h . c . in standard notations, where the nearest-neighbor hopping t i , j = t + t K + δ t i , j is real. We considerthe Kekul´e bond ordering given as t K = ∆ K ( − ∆ K ) for the thick (thin) bonds in Fig. 1 (b), while δ t i , j represents a random component which breaks the translational invariance of the system. Themagnetic field perpendicular to the system is incorporated as the Peierls phase θ i j determined sothat the sum of the phases along a closed loop is equal to the magnetic flux enclosed by the loopin units of the flux quantum ( h / e ). Let us denote the magnetic flux per hexagon as φ .
3. Robust n = Based on the tight-binding model on the honeycomb lattice, we examine the robustness ofthe n = n = n = δ t i , j that obeys a Gaussian probability distribution with avariance σ , P ( δ t i , j ) = √ πσ exp − δ t σ ! . e e (a) Honeycomb Lattce (b) Kekulé structure y x e’e’ y k x k y (c) Brillouin zones Γ K K’
Figure 1: (Color online) (a) The usual honeycomb lattice, for which a unit cell is indicated by the primitive vectors e = ( √ / , / a and e = ( − √ / , / a with the distance a between the nearest-neighbor sites. The filled (open)circles represent the A (B) sub-lattice. (b) The honeycomb lattice with the Kekul´e bond ordering. The thick (thin) bondrepresents the hopping amplitude t + ∆ K ( t − ∆ K ). The unit cell is now given by the primitive vectors e ′ = (3 √ / , / a and e ′ = ( − √ / , / a . (c) The first Brillouin zone for the honeycomb lattice (the outer hexagon) is folded (arrows)onto the inner hexagon (blue) in the presence of the Kekul´e bond ordering, where the Brillouin zone becomes 1 / Γ = (0 , The random components are assumed to be spatially correlated as h δ t i , j δ t k , l i = h δ t i exp( − ( r i , j − r k , l ) / η ) , where η denotes the spatial correlation length and r i , j denotes the spatial position of the bond t i , j .We should note that the present disordered system still respects the chiral symmetry.We evaluate with the Green function method [15] the averaged density of states ρ ( E ) = −h Im G ii ( E + i ε ) i /π with G ii ( E ) = h i | ( E − H ) − | i i , where the angular bracket hi denotes theaverage over the sites i . If we look at the result in Fig. 2. the split n = E = ± ∆ K exhibit an anomalous sharpness as soon as the spatial correlation length η exceeds thenearest-neighbor distance a of the honeycomb lattice. Thus we can conclude that the anomalousrobustness of the n = n = ff ectivelow-energy Hamiltonian with the in-plane Kekul´e distortion, which is given by [17] in a basis(K A , K B , K ′ A , K ′ B ) as H e ff = γ ( k x σ x + k y σ y ) − i ∆ σ y i ∆ σ y γ ( k x σ x − k y σ y ) ! = γ ( k x − ik y ) 0 − ∆ γ ( k x + ik y ) 0 ∆ ∆ γ ( k x + ik y ) − ∆ γ ( k x − ik y ) 0 , where σ ’s are Pauli matrices, ∆ ≡ ∆ K represents the Kekul´e-induced mixing between the K andK’ points, and γ = at / B , the e ff ective Hamiltonianbecomes H e ff = γπ − ∆ γπ † ∆ ∆ γπ † − ∆ γπ . -0.5 h/a D e n s i t y o f S t a t e s E/t
Figure 2: (Color online) Density of states ρ ( E ) = −h Im G ii ( E + i ε ) i /π for the honeycomb lattice with the Kekul´e bondordering ∆ K / t = .
05 in a uniform magnetic field φ/ ( h / e ) = /
50 is plotted for various values of the spatial correlationlength η for the random component of the hopping amplitudes. We assume the strength of disorder to be σ/ t = . ε/ t = . × − . Here π = ( p x − ip y ) / ~ , where p denotes the dynamical momentum p = − i ~ ∂ + e A , and A standsfor the vector potential with ∇ × A = B . The Hamiltonian remainss chiral-symmetric, since Γ H e ff Γ = − H e ff with Γ = σ z σ z ! . For ∆ =
0, the n = H e ff are exact zero-energy states as given by ψ K0 = φ and ψ K ′ = φ with πφ =
0. Note that any linear combinations of those zero-energy eigenstates again belongto the n = ∆ , H ψ ± = ± ∆ ψ ± . where ψ + = φ φ and ψ − = φ − φ πφ =
0. As long as the bond disorder is spatially correlated over several lattice constantsand hence has little e ff ect on the valley mixing, its e ff ect must appear in the gauge degrees offreedom alone. In such a case, the argument due to Aharonov and Casher [23] may be applied inthe same way as in the case of ∆ =
0, which implies that those Landau levels are not broadenedby disorder.
4. Zero-modes around topological defects
Next, we consider a topological defect (vortex) of the Kekul´e bond ordering on the honey-comb lattice, where the topological zero-modes around the vortex emerge in the bulk energy gap[8, 9, 11, 12] irrespective of the presence or absence of uniform magnetic fields perpendicular tothe system. In Ref.[19], topological interface states localized along domain boundaries betweendi ff erent realizations of the Kekul´e bond ordering have been discussed. If we have three or morerealizations, we can have vortices as in Fig.3. For vortices, the number of topological modes at E = lattice sites. Inthe implementation, we take the Chebyshev polynomials up to 8200-th order with the Jacksonkernel for the evaluation of the local density of states [16].Taking advantage of the large system sizes, we examine the robustness against disorder ofthe topological zero-modes not only for a single vortex with the winding number n w = n w >
1. We also examine the e ff ect of disorderfor the topological zero modes for a pair of a vortex and an anti-vortex. In actual numericalcalculations, we place the topological structure in the central region of our system and adopt theopen boundary condition at the boundary. As we shall see in the following, the boundary e ff ectis negligible as long as we discuss the topological zero-modes localized in the central region ofthe system. We first study the robustness of the zero-mode generated by a single vortex with the simplestwinding number n w = ff ect of spatial correlation of the bonddisorder δ t i , j , we again assume that δ t i , j is Gaussian-distributed with a variance σ with the spatialcorrelation length η as in the previous section. The numerical results with the kernel polynomialmethod are shown in Figs 3(b), 4 and 5. Figure 3 (b) shows that the zero-mode wave functionis spatially localized around the vortex, and resides only on the A sub-lattice. The local densityof states at the center site of the vortex in a magnetic field is shown in Fig. 4. It is clearly seenthat, irrespective of whether the spatial correlation of the bond disorder is long-ranged ( η/ a = η/ a = E =
0, suggesting that the energy of the zero-mode is una ff ected by the bond disorder. The resultalso suggests that the spatial correlation of disorder is insignificant for the robustness of thezero-modes induced by the vortex. We therefore confine ourselves to the spatially uncorrelateddisorder for the fermion-vortex system in the following. To confirm the accuracy of the presentnumerical approach with the kernel polynomial method, we also evaluate the local density ofstates well away from the vortex center, for example in the red shaded region in Fig. 3 (a). Insuch a region, the e ff ect of the vortex is hardly seen and the Landau levels for the uniform Kekul´ebond ordering are recovered [Fig. 5] . In particular, the anomalous sharpness of the n = a) (b) Figure 3: (Color online) (a) A vortex structure of the Kekul´e bond ordering with the winding number n w =
1. Hoppingamplitudes with t + ∆ K are indicated by thick blue, green and red bonds, while those with t − ∆ K by thin black bonds.The Kekul´e patterns in the red, green, and blue regions are shifted from each other by a single primitive vector of theoriginal honeycomb lattice. The mismatch of the three Kekul´e patterns creates a topological defect, a vortex. The vortexis assumed to be located at the center of the system which has 10 sites in the present study. (b) The amplitudes ofthe zero-mode wave function ψ ( r ) around the vortex, where the radius of each purple circle represents | ψ ( r ) | and thezero-mode resides only on the A sub-lattice for this vortex. The | ψ ( r ) | satisfies P r | ψ ( r ) | = ∆ / t = . levels is seen to be retained only for long-ranged disorder ( η/ a =
2) in Fig. 5(b). This suggeststhat the present approach has su ffi cient accuracy for discussing the robustness of the zero-modes. Let us move on to a pair of a vortex with the winding number n w = + n w = − ∆ E decays exponentially withthe distance d between the vortices as ∆ E ∝ exp( − d /ξ ). The decay length is estimated to be ξ ≃ a for the case of ∆ / t = .
15, which is in good agreement with the localization length ofthe zero-mode wave function given as γ/ ∆ [8] for the e ff ective Hamiltonian H e ff .We then examine the e ff ect of disorder by adding a spatially uncorrelated bond disorder. Thelocal density of states at the center site of the vortex is shown in Fig. 7 for various values of thedistance between the vortices, d . We can clearly see that the disorder degrades the sharpness ofthe split zero-modes when the vortices are close to each other. The sharpness is recovered as d is increased, where the splitting becomes small. The result indicates that the distance betweenvortex and anti-vortex has to be larger than the localization length of the zero-mode for therobustness of the zero-mode. So far we have considered the simplest winding number n w =
1. How about the vorticeshaving higher winding numbers? So let us look at the result with the winding number for n w = a) (b) η/ a =0 η/ a =2 Figure 4: (Color online) The local density of states ρ ( E ) at the center of the vortex with ∆ / t = .
15 depicted inFig. 3, evaluated by the kernel polynomial method. The uniform external magnetic field perpendicular to the systemis φ/ ( h / e ) = /
50, and the strength of the bond disorder is σ/ t = .
1. An average over 100 realizations of disorder isperformed. Results for (a) η = η/ a = (a) (b) η/ a =0 η/ a =2 Figure 5: (Color online) The local density of states ρ ( E ) ( ∼
200 sites) away form the center of the vortex in Fig. 3, withthe parameters the same as in Fig.4. -1.0 -0.5 0.0 0.5 1.00123 L o ca l D e n s it yo f S t a t e s E / t Δ E (b)(a) Figure 6: (Color online) (a) A vortex with n w = n w = − d in real space. (b) The local density of states ρ ( E ) at the center site (A sub-lattice) of the left vortex forthe case of d / a = ∆ / t = .
15, and the magnetic fluxper hexagon is φ/ ( h / e ) = / ∆ E as a function of the distance d between thevortices. d / a = 4 d / a = 10 d / a = 28 Figure 7: In the presence of a bond disorder the local density of states ρ ( E ) evaluated at the center site of the vortexillustrated in Fig.6 (a) is shown for the inter-vortex distance d / a =
4, 10 and 28. The disorder is here assumed to bespatially uncorrelated with a uniform distribution over [ − W , W ] with W / t = .
2. Other parameters are the same as in Fig.6. An average over 50 realizations of disorder is preformed. P w = 2 w = 4(a) (b) Figure 8: The local density of states ρ ( E ) at the site P ( ∈ A sub-lattice) located in the central region of the vortex with (a)the winding number n w =
2, and (b) n w = ∆ / t = .
6, and a bond disorder, spatially uncorrelated, is considered with a uniformdistribution over [ − W , W ] with W / t = .
2. Here a result for a single realization of disorder is presented. [Fig. 8(a)] and n w = ff ect of n w . The localdensity of states in the central region of the vortex evaluated by the kernel polynomial methodis shown in the figure. We can clearly see that the energy of the zero-modes is una ff ected bythe bond disorder for these higher winding numbers. By evaluating he total amplitude aroundthe vortex, we confirm that the number of zero modes at E = E = ff erent bond patterns in the central region of the vortex. We then find thatthe zero-modes are insensitive to the detailed structure of the central region of the vortex, whichsuggests that the emergence of the zero-modes at the vortex is indeed dictated by the non-localtopological property, such as the winding number, of the vortex structure. All the zero-modesagain reside only on the A sub-lattice, and are still the eigenstates of the chiral operator with theeigenvalue 1. What will happen when the chiral symmetry is broken, typically by a staggered potentialdefined as the Hamiltonian, H µ = µ X i ∈ A c † i c i − X i ∈ B c † i c i . When the chiral symmetry is degraded, irrational charges associated with vortices have beendiscussed [8, 9, 11]. Here we can examine whether the irrational charges remain robust againstthe bond disorder with direct numerical calculations. Even in the presence of the bond disorderwhich respects the chiral symmetry, the energy of the zero-mode associated with the vortexmust be exactly E = µ , since the zero-mode has amplitudes only on the A sub-lattice, hencean eigenstate of H µ . This is a direct consequence of the chiral symmetry for µ = W/t = 0.0
W/t = 0.2
W/t = 0.4
W/t = 0.6 Q / e m/t Figure 9: The fractional charge, Q , associated with the vortex with winging number n w = µ for various values of the strength of disorder W . Spatially uncorrelated bond disorder with a uniformdistribution in the range [ − W / , W /
2] is adopted. The curve represents the formula ( n w / − µ/ p µ + ∆ ) in Ref. [9]. The charge of a vortex is defined as the increment in the local charge when the vortex is in-troduced [8]. We evaluate the charge for various values of the disorder strength and the staggeredpotential µ for the vortices with w = −
4. The result indicate that the charge is proportional tothe winding number n w , and insensitive to the strength of the bond disorder. This is exemplifiedfor n w = W in Fig. 9. The insensitivity of the irrational chargeto disorder may be related to the fact that the bond disorder does not change the energy E = µ ofthe zero modes.
5. Summary
We have investigated numerically the honeycomb lattice with the Kekul´e bond ordering to ex-amine the robustness of zero-modes against disorder that respects the chiral symmetry. We haveclearly shown the following: (i) The split n = n = ff ective theory [9]. The fractional charge is insensitive to the bond disorder that espects thechiral symmetry. It is an interesting future problem how we can relate these with experimentallyobservable properties. 10 cknowledgements This work is partly supported by JSPS KAKENHI Grants Numbers JP19K03660 (TK), andJP17H06138.
References [1] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, A. A.Firsov, Nature , 197 (2005).[2] Y. Zhang, Y-W. Tan, H. L. Stormer, P. Kim, Nature , 201 (2005).[3] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, A. K. Geim, Rev. Mod. Phys. , 109 (2009).[4] S. Das Sarma, S. Adam, E. H. Hwang, E. Rossi, Rev. Mod. Phys. , 407 (2011).[5] V. N. Kotov, B. Uchoa, V. M. Pereira, F. Guinea, A. H. Castro Neto, Rev. Mod. Phys. , 1067 (2012).[6] T. Ando, J. Phys. Soc. Jpn. , 777 (2005).[7] Physics of Graphene , edited by H. Aoki and M. S. Dresselhaus (Springer, Heidelberg, 2014).[8] C.-Y. Hou, C. Chamon, and C. Mudry, Phys. Rev. Lett. , 186809 (2007); Phys. Rev. B , 075427 (2010).[9] C. Chamon, C.-Y. Hou, R. Jackiw, C. Mudry, S.-Y. Pi, and A. P. Schnyder, Phys. Rev. Lett. , 110405 (2008).[10] C. Chamon, C.-Y. Hou, R. Jackiw, C. Mudry, S.-Y. Pi, and G. Semeno ff , Phys. Rev. B , 235431 (2008).[11] S. Ryu, C. Mudry, C.-Y. Hou and C. Chamon, Phys. Rev. B , 205319 (2009).[12] R. Jackiw and P. Rossi, Nucl. Phys. B , 681 (1981).[13] K.K. Gomes, et al., Nature , 306 (2012).[14] Y. Yang, Y. F. Xu, T. Xu, H.-X. Wang, J.-H. Jiang, X. Hu, and Z. H. Hang, Phys. Rev. Lett. , 217401 (2018).[15] L. Schweitzer, B. Kramer, and A. MacKinnon, J. Phys. C , 4111 (1984).[16] A. Weisse, et al., Rev. Mod. Phys. , 275 (2006).[17] N. A. Viet, H. Ajiki and T. Ando, J. Phys. Soc. Jpn. , 3036 (1994).[18] H. Ajiki and T. Ando, J. Phys. Soc. Jpn. , 260 (1994).[19] Y. Hatsugai, T. Fukui, and H. Aoki, Physica E , 1530 (2008).[20] K. Nomura, S. Ryu, and D.-H. Lee, Phys. Rev. Lett. , 216801 (2009).[21] T. Kawarabayashi, Y. Hatsugai, and H. Aoki, Phys. Rev. Lett. , 156804 (2009); Physica E , 759 (2010); Phys.Rev. B , 165410 (2012).[22] T. Kawarabayashi, T. Morimoto, Y. Hatsugai, and H. Aoki, Phys. Rev. B , 195426 (2010).[23] Y. Aharonov and A. Casher, Phys. Rev. A , 2461 (1979)., 2461 (1979).