Berry Curvature and Nonlocal Transport Characteristics of Antidot Graphene

Antidot graphene denotes a monolayer of graphene structured by a periodic array of holes. Its energy dispersion is known to display a gap at the Dirac point. However, since the degeneracy between the A and B sites is preserved, antidot graphene cannot be described by the 2D massive Dirac equation, which is suitable for systems with an inherent A=B asymmetry. From inversion and time-reversal-symmetry considerations, antidotgraphene shouldtherefore have zeroBerrycurvature.Inthiswork,we derive theeffective Hamiltonian of antidot graphene from its tight-binding wave functions. The resulting Hamiltonian is a 4 × 4 matrix with a nonzero intervalley scattering term, which is responsible for the gap at the Dirac point. Furthermore, nonzero Berry curvature is obtained from the effective Hamiltonian, owing to the double degeneracy of the eigenfunctions. The topological manifestation is shown to be robust against randomness perturbations. Since the Berry curvature isexpectedto induce a transverse conductance, we have experimentally verified this feature through nonlocal transport measurements, by fabricating three antidot graphene samples with a triangular array of holes, a fixed periodicity of 150 nm, and hole diameters of 100, 80, and 60 nm. All three samples display topological nonlocal conductance, with excellent agreement with the theory predictions.

To clarify the topological properties of antidot graphene, the first step is to obtain its effective Hamiltonian.Usually, the effective Hamiltonian is ℏν F σ • k þ mσ z [36][37][38][39], which gives the hyperbolic dispersion relation that can fit the band structure around the Dirac point well.Such a Hamiltonian is a 2 × 2 matrix, indicating that there is no valley mixing.The mσ z term models the potential difference between the two sublattice sites [39,40].Since there is no such sublattice asymmetry in antidot graphene, it follows that the massive 2D Dirac Hamiltonian ℏν F σ • k þ mσ z cannot be the effective Hamiltonian for antidot graphene.In the present work, we exploited the tight-binding wave functions and band structures to derive the effective Hamiltonian for antidot graphene in the vicinity of the Dirac point.The new Hamiltonian turns out to be a 4 × 4 matrix that preserves the sublattice symmetry as expected.A similar m term, responsible for opening a gap, appears, which describes the intervalley scattering strength.Since the eigenfunctions are doubly degenerate, nonzero Berry curvature can exist in this system while preserving both time-reversal and inversion symmetries.
A nonzero Berry curvature is expected to induce a transverse conductance [5,9,41,42].We have experimentally verified the existence of this nonzero Berry curvature by measuring the nonlocal transverse current [6,11] and compared the results with the theory prediction.Excellent agreement is found.
A very important point about the opening of a gap in the dispersion relation and the accompanying topological properties of antidot graphene is that they are robust against randomness perturbations.This accounts for their experimental observations in plasma-etched samples.
In what follows, derivation of the effective antidot graphene Hamiltonian is described in Sec.II, followed by the evaluation of the Berry curvature in Sec.III.The measurements of nonlocal transverse electrical resistance arising from the Berry curvature, attendant with the simulation of the theory prediction, are presented in Sec.IV.Comparison between the measured data on three fabricated samples and the theory predictions is given in Sec.V. Section VI concludes with a brief summary of results.

II. EFFECTIVE HAMILTONIAN OF ANTIDOT GRAPHENE A. Tight-binding model
Consider a triangular lattice of holes (antidots) in monolayer graphene with lattice constant L as shown in Fig. 1(a) .The unit cell is hexagonal in shape with a hole at its center.To model antidot graphene, we define a dimensionless geometrical factor γ ¼ d=L, where d denotes the hole diameter.Therefore, (L, γ) can be used to describe the antidot unit cell.The band structure and wave function can be calculated by using the tight-binding model, where the Hamiltonian is given by Here a † i denotes the creation operator on the ith A sublattice site, b iþδ denotes the destruction operator on the i þ δth B sublattice site, and δ includes the three nearest neighbors to site i.In Eq. (1), the summation over i goes over the whole unit cell.We apply the Bloch boundary condition at the edge of the hexagons.The band structure and wave function of the antidot graphene system can be obtained by directly diagonalizing the Hamiltonian represented by Eq. ( 1).The obtained band structure is shown by the black circles in Fig. 1(c) for various γ with a fixed periodicity L ¼ 13 nm.The band structure is plotted as a function of the Bloch wave vector k, which ranges from K to Γ, then to M, and finally back to K. Here, the K, Γ, and M points are for the hole reciprocal lattices.As the antidots are arranged in the triangular lattice structure, its first Brillouin zone is hexagonal in shape, the same as the first Brillouin zone of the pristine graphene lattice.The Dirac point, i.e., the K point of the pristine graphene's reciprocal lattice, folds onto the Γ point of the hole reciprocal lattice.We can see that a gap opens at the Dirac point.Both the conduction and valence bands are doubly degenerate because of the existence of K and K 0 valleys.We have chosen to ignore some flat bands [43] since they correspond to localized states and do not contribute to the transport properties.At the same time, we can also obtain the wave functions ψ ¼ ðφ A ; φ B Þ T , which can be analyzed in the k domain by carrying out the Fourier transform.

B. Effective Hamiltonian
Since we wish to focus on the low-energy region, the Fourier coefficients for both φ A and φ B are largely determined by two components, i.e., around the K and K 0 .In other words, and similarly for the wave function on the B site, Therefore, the wave function can be represented by a four-dimensional vector ða K ; a K 0 ; b K ; b K 0 Þ T , with bases e iðKþkÞ•r and e iðK 0 þkÞ•r for both A and B atoms, respectively.It follows that if we define where i ¼ 1 or 2 denotes the two conduction-band wave functions, and i ¼ 3 or 4 denotes the two valence-band wave functions; then the effective Hamiltonian H should satisfy the following equation: PAN, ZHANG, ZHANG, ZHANG, DONG, and SHENG PHYS.REV.X 7, 031043 (2017) where E þ and E − denote the (degenerate) eigenvalues of the corresponding eigenfunctions and I is the 2 × 2 identity matrix.The effective Hamiltonian can be obtained by With the eigenvalues and eigenfunctions obtained from the tight-binding calculation, we have numerically evaluated the effective Hamiltonian, which is in the form of a 4 × 4 matrix with only eight nonzero matrix elements: We have found that jH 12 j ¼ jH 21 j ¼ jH 34 j ¼ jH 43 j and jH 14 j ¼ jH 23 j ¼ jH 32 j ¼ jH 41 j.In Fig. 1(b), we present the results of the matrix element as a function of k around the Γ point of the antidot lattice.It is seen that jH 12 j is linear in k, while jH 14 j is almost a constant.The details of H 12 's k dependence (i.e., as a function of k x and k y ) can be obtained by setting k along the x and y directions.It is found that the Hamiltonian can be expressed as Therefore, the effective Hamiltonian consists of two 2 × 2 matrices, ℏvσ • p and −ℏvσ Ã • p, along the diagonal and four antidiagonal m terms [44].The two diagonal matrices are seen to originate from the pristine graphene Hamiltonian, which can be easily verified by doing Taylor expansion around K and K 0 points of pristine graphene.Note that the velocity v here is not necessarily the same as v F (c=300) in pristine graphene.In antidot graphene, it is a function of the geometric factor γ. The four antidiagonal m terms couple the K and K 0 valleys; hence, they imply intervalley scatterings, which can arise from the atomically sharp edges of the antidots.We find the intervalley scattering strength to have a very weak k 2 dependence, which is ignored here since we focus only on the low-energy regime.Hence, the intervalley scattering strength is characterized by a constant m.The effective Hamiltonian yields a dispersion relation , where both the conduction and valence bands are doubly degenerate.This dispersion relation, which is shown as red curves in Fig. 1(c), can fit the numerically obtained band structures very well, thereby confirming this effective Hamiltonian to be relatively accurate for modeling the antidot graphene system in the vicinity of the Γ point (of the antidot lattice).
The values of the velocity, v and m, depend on the geometry of antidot graphene.The values evaluated from the tight-binding calculations are summarized in Fig. 2. For a fixed periodicity L, the effective velocity v (black solid squares) decreases with increasing γ, while the gap (2m) increases with increasing γ.A decrease of the Fermi velocity associated with the gap opening can also be found in Ref. [45].It can be plausibly understood as due to the narrowing of the "neck" region of the passage channels with increasing γ.We can linearly fit m and v as a function of γ, which are shown as dashed curves in Fig. 2: The term m scales with the inverse of the periodicity, i.e., varying as 1=L, which is in accordance with previous numerical results [19,46].For our experimental samples, we have a fixed L ¼ 150 nm and γ ¼ 2=3 (sample A), 8=15 (sample B), and 2=5 (sample C).The corresponding gap (2m) is around 20 meV, and the effective velocity v is about 0.55v F ; hence, m=ℏv ¼ 0.03 nm −1 for sample A. Values of m and v for samples B and C can be similarly deduced from Eq. (8).

C. Uniqueness of the effective Hamiltonian
Since the bands are doubly degenerate, the two wave functions ð ψ 1 ψ 2 Þ at the same energy can form a linear combination ð ψ 0 1 ψ 0 2 Þ, which would be an equally valid choice.The question is whether different choices of the wave functions can result in a different effective Hamiltonian.To prove the uniqueness of the effective Hamiltonian, we recognize the fact that the wave functions must transform into each other via the unitary matrix U AE for both the conduction and valence bands.In other words, Let H 0 be the alternative effective Hamiltonian and let Based on Eq. ( 5), we have Here, Therefore, we can conclude that the effective Hamiltonian obtained above is unique; i.e., it would remain the same regardless of the choice of wave functions.

A. Symmetry analysis
The nonzero term m mixes the K and K 0 valleys whose Berry curvatures have opposite signs.Hence, one may wonder whether there is net nonzero Berry curvature with this valley mixing.It is noted that both time-reversal and inversion symmetries are preserved in this system; however, that does not guarantee a zero Berry curvature [9,41,42,47,48].We start the symmetry analysis with the eigenvectors of the conduction band, obtained directly from the effective Hamiltonian, Eq. (7),

<
: where θ ¼ arctanðk y =k x Þ denotes the direction of the wave vector k.Since there is double degeneracy for the conduction band, any linear combination of ψ 1 and ψ 2 is also an eigenwavefunction of the effective Hamiltonian, Eq. (7).In Eq. ( 12), the eigenwavefunctions are chosen so that the Berry curvature matrix is diagonalized, as shown in Sec.III B. Under time-reversal transformation (see Appendix A), Tψ 2 ðkÞ ¼ ψ 2 ð−kÞ Ã ¼ e iθ ψ 1 ðkÞ; ð13Þ from which we conclude that the Berry curvature must obey (see Appendix B) where the subscript denotes the respective wave function [Eq.( 12)] from which the Berry curvature is evaluated.
Similarly, for the inversion symmetry, we have Ω 1 ð−kÞ ¼ Ω 1 ðkÞ and Ω 2 ð−kÞ ¼ Ω 2 ðkÞ.Combining both the timereversal and inversion symmetries, we find that the symmetry conditions only require Hence, the Berry curvature for each band is not necessarily zero.Below, we can see that, for our degenerate system, the Berry curvature is indeed nonzero.
If we focus on the region near the Γ point (of the hole lattices), i.e., jkj ∼ 0, the eigenvectors from Eq. ( 12) can be simplified as Equation ( 16) clearly shows the mixing behavior between the K and K 0 valleys.For instance, for ψ 1 , the wave function of atom A comes from the K 0 valley, while for atom B, it comes from the K valley.However, if we assume that there is a long-range potential V (such as the Coulomb potential of the charged impurity), the scattering matrix element should be In other words, the long-range disorders cannot scatter ψ 1 to ψ 2 .Note that the intervalley scattering is already taken into account in the Hamiltonian and its relevant wave functions.

B. Calculation of the Berry curvature
Since the 4 × 4 matrix describes a doubly degenerate system, the Berry curvature for either the conduction or valence band is essentially a 2 × 2 matrix, which is defined as [9,41,42,[47][48][49]] Here, the Berry connection A α is a matrix, whose matrix elements are defined as ðA α Þ mn ¼ ihu m j∂k α ju n i, where u mðnÞ denotes the periodic part of the wave function.We can see that the first two terms originate from the nondegenerate case, while the last term is the correction term due to degeneracy [9,42,[47][48][49].If there is no degeneracy, A α is a number instead of a matrix, making ½A x ; A y zero.
The calculation of the Berry connection requires information from the periodic parts of the wave function denoted by u mðnÞ .We can derive the expression for u based on Eqs.(2a) and (2b , where ψ mðnÞ are the eigenvectors of the effective Hamiltonian [Eq.( 12)].Hence, by substituting Eq. ( 12) into Eq.( 18), we find that the Berry curvature is given by It is seen that the Berry curvature matrix is nonzero only for the diagonal terms.It is also noted that the matrix element Ω 11 is the same as the Berry curvature of the massive Dirac Hamiltonian [5], where m characterizes the local potential difference between the carbon atoms A and B. The Berry curvature is determined solely by one parameter, i.e., m=ℏv.In Fig. 3, we summarize the Berry curvature for different geometric parameters γ.It peaks at k ¼ 0 and decays as k deviates from 0. The Berry curvature can also be obtained through tight-binding calculations [48].We found that near the Γ point of the antidot lattice, Eq. ( 19) agrees well with the Berry curvature obtained through tight-binding calculations.Away from the Γ point, there can be differences between the two.It is also to be noted that by integrating the Berry curvature calculated by the tight-binding model over the first Brillouin zone, we obtain a Berry phase of 2π.Therefore, the effective Hamiltonian, Eq. ( 7), and the Berry curvature, Eq. ( 19), can well describe the antidot graphene system only at the low-Fermi-energy region.

C. Robustness against randomness perturbations
From the derived effective Hamiltonian and the related analysis, it should be clear that the gap opening in the antidot graphene is due to the folding of the K and K 0 points [26] of the reciprocal carbon lattice, onto the Γ point of the hole reciprocal lattice.Specifically, K and K 0 (of carbon lattice) are located at (AE4π=ð3 ffiffi ffi 3 p aÞ; 0), where a is the carboncarbon atomic separation.If the periodicity of holes is given by 3Na, where N is an integer, then the primitive vector of the reciprocal hole lattice is 4π=ð3 ffiffi ffi 3 p NaÞ • ð1=2; AE ffiffi ffi 3 p =2Þ.We can see that K and K 0 can be folded onto the Γ point.
In the above context of a perfect antidot lattice, it was pointed out in Ref. [22] that if the hole lattice is rotated by π=6, which is denoted the "rotated triangular lattice," and, in addition, if the periodicity of the hole lattice is ffiffi ffi 3 p Na so that the primitive vector of the reciprocal hole lattice is 4π=ð3NaÞ • ð ffiffi ffi 3 p =2; AE1=2Þ, with N ≠ 3n (n being an integer), then K and K 0 valleys are not coupled; therefore, there can be no band gap.These conclusions were numerically verified in Ref. [22], where the results were explained by using the Clar sextet theory.Below, we show that such cases are unstable against small-scale randomness perturbations on the edge of the holes.Such randomness would be difficult to avoid in the oxygen etching of the antidot pattern.It is demonstrated that once there is atomic-scale randomness on the edges of the holes, then the full band gap is restored; the magnitude of the band gap in such cases is fairly accurately predicted by Eq. ( 8) from the geometric parameter γ.
The conclusion from such randomness perturbation analysis, shown below, is that the results obtain from Sec. II are generally applicable to the experimentally fabricated antidot samples.
To determine whether the antidot sample is gapped or gapless, we calculate its conductance as a function of Fermi energy.If there is a zero conductance plateau, then there is a band gap; otherwise, it is a gapless system.The conductance of the antidot graphene system can be well simulated by the open-source package KWANT [50], which can calculate the electronic transport property based on the tight-binding models.We have carried out simulations for the three different cases as follows.(a) Lattice A-This is the "triangular lattice structure," which is what we have focused on in Sec.II.In this case, K and K 0 can always be folded onto the Γ point of the hole reciprocal lattice; i.e., a gap is expected.
Here, the periodicity is fixed at L ¼ 24a, with γ ¼ 0.2.(b) Lattice B-This is the "rotated triangular lattice structure" [22], which is predicted to have no gap.The periodicity in this case is fixed at L ¼ 14 ffiffi ffi 3 p a with the same γ ¼ 0.2 as in (a).The value of the periodicity, L ¼ 14 ffiffi ffi 3 p að≈24.2aÞ, is chosen so that K and K 0 cannot be folded onto the same Γ point.A gapless system is expected; i.e., the conductance should exhibit no zero-conductance region.(c) Lattice C-This is the same as Lattice B but modified with atomic-scale randomness at the edges of the holes.We assume that the hole radius is angular dependent; i.e., r n ðθÞ ¼ r where n refers to the nth hole.The term f n ðθÞ is defined as where the summation of m goes from 1 to 5. In the simulations, we have randomly chosen a nm ∈ ð0; 0.06Þ, φ nm ∈ ð0; 2πÞ for each hole.The whole system for transport simulations is shown in Fig. 4(a).We can see that there are 10 holes in the width  A, where the zero-conductance plateau indicates the band gap of this system.For Lattice B, the red curve exhibits no zeroconductance regions, implying that there is no band gap.Lattice C differs from Lattice B only at the edges of the holes, where we have added atomic-scale randomness.Ten different randomness profiles were calculated, and the blue curve represents the ensembleaveraged conductance.It is seen that a band gap reappears, and the magnitude of this band gap is the same as that of Lattice A, implying that the transport gap in Lattice C has the same origin as that in Lattice A, i.e., from K to K 0 intervalley scatterings.direction and 20 holes in the length direction.The inset image is a zoom-in for a hole in Lattice B, i.e., "the rotated triangular lattice structure."The conductances as a function of Fermi energy of the three cases are summarized in Fig. 4(b).For Lattice A, the conductance is plotted as the black curve, where a zero conductance plateau is observed as expected.For Lattice B, the conductance does not drop to zero at zero Fermi energy, shown by the red curve.This clearly implies that there is no band gap for Lattice B, as predicted.For the third case, we introduced atomic-scale randomness into Lattice B. In this case, we have calculated 10 different randomness profiles, and the blue curve in Fig. 4(b) is the averaged conductance.It is clearly seen that with randomness, the full band gap reappears.The magnitude of the band gap in the third case is exactly the same as that of Lattice A, around 0.08t.With the appropriate parameter values of L ¼ 24a and γ ¼ 0.2, Eq. ( 8) predicts a band gap around 0.084t, which is consistent with the conductance simulation results.
There is a simple explanation for the reappearance of the gap when the atomic-scale perturbations were added.In this case, small-scale randomness in real space implies scatterers with large momentum transfers in the reciprocal space, i.e., short-range scatterers.That means intervalley scatterings are the inevitable consequences of introducing atomic-scale randomness, and the opening of the gap is the result.There is also a simple heuristic explanation for why the obtained gap size depends on the geometric parameter γ=L ¼ d=L 2 as given by Eq. ( 8).The K to K 0 scattering strength is given by hKjV sr jK 0 i, where V sr denotes the effective short-range disorder potential.Here, V sr must be proportional to the number of sites at the edges of the holes, which is proportional to the circumference of the hole; hence, V sr ∼ nd, where n is the number of holes.The wave amplitudes of jKi and jK 0 i must each be normalized in the sample area S, which means jKi and jK 0 i must each be proportional to 1= ffiffiffi S p .Combining these two facts, we get hKjV sr jK 0 i ∝ nd=S ∼ d=L 2 , as S ∝ nL 2 .However, such an argument cannot yield the accurate value of the gap, which must be determined numerically.
From this perspective, the existence of a gap in antidot graphene is almost inevitable for the experimentally fabricated antidot samples by using plasma etching.The nogap state is very precarious since its existence requires the very precise periodicity, with no atomic-scale randomness allowed.

D. Energy splitting of the eigenfunctions
Accompanying the introduction of edge randomness into the antidot graphene is the splitting of the originally doubly degenerate eigenfunctions as in the case of Lattice A above.Below, we first quantify such splittings and then show in the subsequent section that the topological manifestation of antidot graphene is not affected by such splittings.
We have calculated the band structure of antidot graphene with randomness by using the supercell approach in which each supercell comprised several unit cells with randomness introduced at the hole edges.The Bloch boundary condition was imposed at the boundary of supercells.The periodicity of the unit cell is set to 19 ffiffi ffi 3 p a so that there is no K to K 0 coupling in the perfect antidot lattice case.For the supercell, we consider nine holes as shown in Fig. 5.The periodicity of the supercell is For the band-structure data presented in Fig. 5, we have removed the flat bands.If fðθÞ ¼ 0, i.e., no randomness, we can see that there is no band gap, as shown in Fig. 5(a).For Fig. 5(b), we consider the same case as in Fig. 5(a) but with atomic-scale randomness.It is noted that there is a tiny energy splitting of the originally doubly degenerate eigenfunctions, accompanying the appearance of the band gap.Hence, the conduction (valence) bands become quasidegenerate.In Fig. 5(b), the size of the gap is around 0.05t, which is close to the prediction by Eq. ( 8), which gives a gap around 0.052t.
In the presence of the randomness perturbation, the effective Hamiltonian can be numerically evaluated to be FIG. 5. Supercells and the simulated band structures.Each supercell comprises nine unit cells, whose hole periodicity is given by 19 ffiffi ffi 3 p a. In diagram (a), we have an antidot lattice without any randomness; i.e., periodicity is accurate down to the atomic scale.We can see that there is no band gap, as predicted.(b) Antidot lattice with randomness.For this case, the band gap reappears, accompanied by a small energy splitting for both the conduction and valence bands.In other words, the states in the conduction and valence bands are quasidegenerate.
This nonzero c turns out to be responsible for the energy splitting of the eigenfunctions in the conduction (valence) band.
In fact, one can easily evaluate the dispersion relation from the above effective Hamiltonian to give EðkÞ , which indicates an energy splitting in the eigenfunctions of both the conduction and valence bands.This splitting turns out to be on the order of 0.1m.Hence, both the conduction and valence bands can be regarded as quasidegenerate.
One may wonder what effect this splitting might have on the topological characteristics of antidot graphene.Below, we show that the topological manifestation is not changed from the case when the eigenfunctions are doubly degenerate.

E. Effect of the Berry curvature under an applied electric field
Berry curvature is well known to affect the electronic transport in the presence of an electric field.This is a physical manifestation of Berry curvature.In the present case, the Berry curvature Ω is a 2 × 2 matrix since the conduction and valence bands are (quasi-) double degenerate.Let ðφ 1 ; φ 2 Þ denote the two quasidegenerate states; an arbitrary initial state ψ is expressible as a linear superposition of ðφ 1 ; φ 2 Þ, i.e., ψ ¼ η 1 φ 1 þ η 2 φ 2 , where η 1;2 are the coefficients satisfying jη 1 j 2 þ jη 2 j 2 ¼ 1.It would be convenient to express ψ in a vector form η ¼ ð η 1 η 2 Þ T .In the presence of an applied electric field, the electronic velocity transverse to the electric field direction is given by [47] where E y is the external electric field.Based on the tightbinding model [48], we have numerically calculated the Ω matrix by using the supercell with randomness.It turns out that the off-diagonal terms of Ω are nonzero.For simplicity in evaluating the effect of the Berry curvature under an electric field, we choose to diagonalize Ω by changing to a new basis.If the diagonalization matrix is U, we have . It follows that Eq. (21a) can be simplified as The values of Ω 11 and Ω 22 in the first Brillouin zone are summarized in the inset of Fig. 6.We can see that Ω 11 and Ω 22 have opposite signs, i.e., Ω 22 ¼ −Ω 11 , and their values are very close to those shown in Eq. ( 19).From Eq. (21b), we have v x ¼ ðeE y =ℏÞ • ðjη 0 1 j 2 − jη 0 2 j 2 ÞΩ 11 .This result means that when jη 0 1 j 2 > 1=2, the electronic wave tends to go left as E y is negative (applied in the minus y direction), as illustrated by the red wave in Fig. 6.If we focus on the states that flow leftwards, the average transverse velocity is given by where the denominator is the normalization factor.The summation can be carried out by integration.By assuming that η 0 1 ¼jη 0 22) can be easily evaluated, and the average velocity turns out to be vx ¼ ðeE y =ℏÞ • ðΩ 11 =2Þ.For the range jη 0 1 j 2 < 1=2, the average transverse velocity is vx ¼ −ðeE y =ℏÞ • ðΩ 11 =2Þ.Hence, the effective Berry curvature for inducing the transverse current is half of the diagonal matrix element, i.e., AEΩ 11 =2.
In the above, it can be seen that independent of whether the states are doubly degenerate or quasidegenerate, the Berry curvature in antidot graphene has the same physical manifestation.The atomic-scale randomness at the edge of the holes, in fact, makes the topological effect more robust, i.e., insensitive to the relative rotation of the hole lattice vs the underlying graphene lattice, as well as to the precise values of the hole lattice constant.This predicted "topological Hall effect" is interesting, and we show experimentally, in the following section, by fabricating and measuring three samples with different geometric parameter values, that the topological current indeed exists in antidot graphene.

A. Sample and experimental setup
To experimentally check the nonzero Berry curvature, we have fabricated the antidot graphene samples and studied their electronic transport properties by exploiting the nonlocal measurements.The single-layer graphene was obtained by mechanical exfoliation and transferred onto a silicon wafer with a SiO 2 thickness of 285 nm.The periodic antidots are patterned by e-beam lithography, followed by oxygen plasma etching.The whole sample was etched into standard Hall bar geometry, shown as the cartoon image in Figs.7(a) and 7(b), where the blue region is the antidot graphene, while the grey parts denote the electrodes (10 nm Ti=60 nm Au).The geometric parameters are set as d 1 ¼ 2 μm, d 2 ¼ 4 μm, and W ¼ 2 μm.For the antidot lattice, we have fabricated three different samples.All samples have periodicity fixed at L ¼ 150 nm, while the hole diameters were designed to be d ¼ 100 (sample A), 80 (sample B), and 60 (sample C) nm.Scanning-electronmicroscope images are shown in the inset of Fig. 7(c) for sample A, and in Fig. 14 for sample B (upper panel) and sample C (bottom panel).All measurements were conducted in PPMS (Quantum Design) after in situ annealing at 390 K for 2 hours.We focus on analyzing the data from sample A. Similar results for samples B and C are summarized in Fig. 14.

B. Nonlocal measurements
The nonzero Berry curvature in antidot graphene can be experimentally verified by applying the nonlocal measurements, and the setup is shown in Fig. 7(a).By applying a current I through one pair of Hall electrodes and recording the voltage drop V nl between another nearby pair of electrodes, the nonlocal resistance is defined as The nonlocal measurement setup differs from the traditional four-probe local measurement, as shown in Fig. 7(b).The traditional four-probe local measurement gives the longitudinal resistivity ρ of the antidot graphene; its dependence on gate voltage for various temperatures is shown in Fig. 7(c).We can see that the charge neutrality point (CNP) is located around 7.5 V.The longitudinal resistivity has a strong temperature dependence, which can be well described by a variable range hopping model (VRH) [29,51].The mobility of antidot graphene is usually low compared to pristine graphene [29,52].The resistivity (conductivity) as a function of gate voltage can be described by ρ [29,[52][53][54], where ρ s comes from the short-range scatterings, μ is the mobility controlled by long-range (chargedimpurity) scatterings, C is the capacitance per unit area of SiO 2 , and n 0 is the residue carrier density due to the electron-hole puddle effect [55][56][57].These parameters can be obtained by fitting the room-temperature data.
The fitting details are shown in Fig. 7(d), which gives mobility μ ¼ 1500 cm 2 =ðVsÞ and residue carrier density

C. Theory on nonlocal transverse conductance and simulations
The nonlocal resistance is determined by the topological current that flows in the sample channel.Since the state has a nonzero Berry curvature, it will result in a transverse conductance σ xy , given by [5,6] where k ¼ jkj, f (k) is the Fermi distribution function, and the Fermi level is determined by the carrier density n, which is related to k F by the relation n ¼ k 2 F =π. Here, the factor 1=2 comes from the evaluation of effective Berry curvature as shown in Sec.III E. Substituting Eq. ( 19) into Eq.(24a), we obtain The nonlocal resistance is strongly related to the transverse conductance and may be expressed as [6,11] where α is a coefficient that is independent of the Fermi energy but depends on the sample geometric parameters, W and d 2 , where W denotes the sample width and d 2 the separation between the electrodes as shown in Fig. 7.
Ultimately, the coefficient α is also related to the decay length ξ of the Berry curvature density, as shown below.
An accurate value of α is important for the comparison between the theory prediction and experimental data.To find the behavior of coefficient α, a differential diffusion equation must be solved.We show below that depending on whether ξ is larger or smaller than the sample width W, different behaviors of α can be derived.
The diffusion equation below was originally formulated to describe spin diffusion, but it can be adapted to describe Berry curvature diffusion [11]: where s is the "Berry curvature" (spin) density.The source term Ξðx; yÞ is determined by the external electric field, with a form given by [11] Ξðx; yÞ ¼ − τ s Iσ xy =Wσ Here I denotes the current, and σ xy , σ stand for the transverse and longitudinal conductivities, respectively; τ s is the relaxation time, which is related to ξ by τ s ∼ ξ 2 .The two delta functions in Eq. ( 26b) indicate that the source term is nonzero only at the two edges (y ¼ AEW=2); hence, Eq. (26b) serves as a boundary condition for sðx; yÞ.The induced Berry curvature current along the x direction is given by Here, the quantity ξ 2 =τ s is the intrinsic diffusion constant; it is treated as a constant of the problem.Based on the above definition of the Berry curvature current, we can obtain the nonlocal resistance, given by Eq. ( 25).In Ref. [11], it is shown that in the range W ≪ ξ, the coefficient α can be analytically expressed as α ¼ W expð−d 2 =ξÞ=2ξ.However, in the regime W ≥ ξ, the coefficient α has to be solved numerically.The case of small ξ (as compared with W) is crucial here, as the decay length must be quite small in antidot graphene since the mobility is low.Below, we show that by fitting the experimental data, the decay length was found to be on the order of 450-550 nm for our three antidot graphene samples.They are consistently smaller than the decay length (2 μm) reported in graphene placed on top of h-BN [6].
Following the work in Ref. [11], we have obtained the coefficient α by numerically solving the differential diffusion equation [Eq.(26a)] with its boundary conditions [Eq.(26b)], by using the commercial software COMSOL [58] package.The Berry curvature density sðx; yÞ can be obtained numerically.One simulated example with ξ ¼ W is shown in Fig. 8(a).The color scale indicates the magnitude of the Berry curvature density; its magnitude is seen to decay to zero as x deviates from x ¼ 0. The Berry curvature current density, defined as J s ðxÞ ¼ −ðξ 2 =τ s Þ∂ x R W=2 −W=2 dy • sðx; yÞ, is summarized in Fig. 8(b) for various values of ξ.Here, J s ðxÞ is seen to decay exponentially with x.Hence, the coefficient α should behave as α ≡ f expð−x=λÞ, where f and λ are dependent on ξ.Their functional dependencies are shown in Fig. 9(a).In the range ξ ≫ W, λ ¼ ξ and f ¼ W=2ξ, our simulations exactly reproduce the analytical expression given previously.However, in the regime ξ ≤ W, λ saturates at W=π when ξ approaches 0. That means when ξ is small, the nonlocal resistance decays exponentially as expð−πx=WÞ, instead of expð−x=ξÞ, as previously predicted by the analytical formula.For the parameter f, its value peaks at W=π and decays quickly as ξ decreases.These behaviors clearly imply that the previous analytical expression for α cannot model the properties of the system in the regime ξ ≤ W. The behavior of the coefficient α, as a function of ξ, is summarized in Fig. 9(b).It is seen that the simulation results deviate from the analytical expression in the small ξ=W range.In the following analysis of the nonlocal resistance, we will use the simulation results in Fig. 9(b).
The exponential decay length has been measured on disordered graphene in Ref. [12], where λ was found to be 300 nm with a sample width of 0.9 μm.The relation between λ and the sample width is noted to agree with what we have predicted here.

D. Simulation of the stray current
In order to isolate the nonzero transverse conductance σ xy , it is necessary to exclude the stray current that can exist even when there is no topological current.Since the voltage-probe electrodes are located on opposite sides of the system [as shown in Fig. 7(a)], a voltage drop between the voltage-probe electrodes is expected once a current flows through the sample.This classical diffusive transport behavior can be understood as the stray current effect, which can be approximated by the van der Pauw relation, i.e., R 0 nl ∼ ρ expð−πL=WÞ.This effect can also be accurately calculated by COMSOL simulations, which is shown in Fig. 10, where the color indicates the potential distribution.To clearly observe the nonlocal voltage difference (which is quite small), we intentionally focus on the range of 0.45 V and 0.55 V.By analyzing the potential drop between the neighboring Hall electrode pair, we have calculated the nonlocal resistance R 0 nl that results from the stray current in our system, where ρ is the longitudinal resistivity of antidot graphene.

E. Nonlocal resistance from the Berry curvature
For the nonlocal measurements, we find that by using the nearest-neighbor electrode pair for applying current and detecting voltages (with the electrode separation of 4 μm), we can always observe a peak at CNP for the nonlocal resistance.However, for electrode pairs of other separations (such as 8 or 12 μm), the detected voltage drop is nothing but noise.In Fig. 11, we present the measured nonlocal resistance R exp nl (solid curves obtained from nearest neighbors) as a function of gate voltage at various temperatures.To rule out the stray current effect, we also plot the calculated R 0 nl (dashed curves) in the same figure.For a clearer comparison, the signal of R 0 nl is amplified by a factor (7 for 10 K and 4 for 20 K) at low temperatures.We can see that the measured nonlocal resistance R exp nl exceeds the stray current effect (R 0 nl ), implying that the measured excess must arise from the topological effect.
We summarize the temperature dependence of the nonlocal resistance in Fig. 12, where the black solid squares represent the peak values (backgate fixed at CNP) of measured nonlocal resistance, and the red solid curve stands for the calculated R 0 nl .Based on Eqs. ( 25) and ( 27), we conclude that the peak values follow the relations R nl ∼ ρ 3 and R 0 nl ∼ ρ.Since ρðTÞ near the CNP decreases when the temperature increases [29] [see Fig. 7(c)], R nl decreases faster than R 0 nl as temperature rises.This is clearly shown in Fig. 12, in which the black dashed curve, representing R nl ðTÞ þ R 0 nl ðTÞ, is seen to display a crossover behavior that agrees extremely well with the measured nonlocal resistance of R exp nl .It follows that at high temperatures, R exp nl is dominated by the stray current effect, whereas at low temperatures, the topological nonlocal resistance dominates.For T < 100 K, the measured nonlocal resistance is larger than what is expected from the stray current effect.This means that for T < 100 K, the topological current can induce a detectable nonlocal voltage drop.We define the topological nonlocal resistance as It is shown below that the carrier concentration dependence of this nonlocal resistance can be well predicted by Eq. (28).Apart from the magnitude, there is also a difference between the features of the measured nonlocal resistance R exp nl and the stray current-induced R 0 nl .For the low-temperature data (10 K and 20 K), the measured nonlocal resistance is zero for jV − V CNP j > 2 V as shown in Fig. 11.For the stray current effect (dashed curves), its values decay to zero smoothly.

V. COMPARISON BETWEEN THEORY AND EXPERIMENT
In Fig. 13, we plot the topological nonlocal resistance, R nl ¼ R exp nl − R 0 nl , as a function of carrier density (open circles).The measured data can be well fitted by the following equation, where the term n in Eq. ( 25) is replaced by ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n 2 þ n 2 0 p , with n 0 ¼ 2 × 10 11 =cm −2 being the residue carrier density at CNP.Since m=ℏv ≈ 0.03 nm −1 is known from the tight-binding calculations, the fitting requires just one parameter-the decay length ξ.The best fit gives ξ ≈ 450 nm.Note that this value is indeed much smaller than the sample width W ¼ 2 μm.
The small decay length values in our samples, as compared to that reported in Ref. [6], can be attributed to the fact that our graphene samples are nanostructured and placed on top of SiO 2 , while in Ref. [6] the sample is pristine graphene placed on h-BN substrate.The latter has been shown to be a much cleaner and flatter substrate [59].In our samples, the magnitude of the residue carrier density is on the order of 10 11 =cm 2 , much larger than that for pristine graphene on h-BN, which is on the order of 10 10 =cm 2 .It follows that our samples are more disordered, and therefore, they can have a much shorter scattering time τ s .Since ξ ∝ ffiffiffiffi τ s p , a short decay length follows.
Based on the decay length ξ ≈ 450 nm, it is easy to understand why we cannot detect the nonlocal voltage drop for the next-nearest neighbor electrode pairs (with a separation of 8 μm).Since the magnitude of nonlocal voltage drop is determined by the coefficient αðξ; xÞ, defined as α ≡ fðξÞ expð−πx=WÞ.For the nearest-neighbor electrode pairs, x ¼ 4 μm, which gives expð−πx=WÞ ¼ 1.8 × 10 −3 .However, for the next-nearest-neighbor electrode, x ¼ 8 μm, which yields expð−πx=WÞ ¼ 3.5 × 10 −6 .Consider the data at 10 K (see Fig. 11); for x ¼ 4 μm, the peak of nonlocal resistance is about 10 kΩ.Based on this value, for x ¼ 8 μm, the peak value of the nonlocal resistance is estimated to be on the order of 10 Ω, which is below the noise floor.
The nonlocal resistance as a function of carrier density measured at 10 K for samples B and C are summarized in   A recent experimental work [60] demonstrated that a conducting edge state might exist in gapped bilayer graphene.Can such a conducting edge state exist in our sample?For our antidot graphene samples placed on top of SiO 2 , the evidence seems to indicate that there cannot be a conducting edge-state channel since the existence of such a state should enhance the measured transverse nonlocal resistance, and over a longer distance, e.g., over 8 μm.These were not seen.More conclusively, we have measured the I-V curve of sample C at 10 K with gate voltage fixed at CNP.The I-V curve is shown in Fig. 15.The nonlinear behavior of the I-V curve indicates the existence of a gap in antidot graphene.From the I-V curve, the derived resistance per square as a function of bias voltage is also plotted.Note that the resistance at CNP (measured under a small bias voltage) is around a few hundred kilo Ohms, which is much larger than the edge-state channel where the resistance should be on the order of h=e 2 ∼ 25.9 kiloOhms.Thus, the existence of a conducting edge state is highly unlikely.

VI. CONCLUSIONS
We want to show that This can be achieved in the following two steps.
(1) By recalling that ψ 2 ðkÞ ¼ f , where the four components correspond to KA, KB, K 0 A, K 0 B. By flipping the sign of k, we must remember that we are changing the total momentum of (K þ k) to −ðK þ kÞ.Since −K ¼ K 0 , it follows that changing the sign of k is always accompanied by exchanging the K and K 0 valleys.Hence, KA becomes K 0 A and KB becomes K 0 B, and vice versa.Therefore, we have FIG. 15.The black curves represent the I-V curve at 10 K for sample C. The red curve is the derived resistance per square.The resistance at small bias voltage is around 100 kΩ, which is much larger than h=e 2 .Hence, the existence of a conducting edge state is not likely.
(1) We first show that the additional phase exp (iθ) does not change the Berry curvature.By replacing ψ 1 in Eq. (B1) with ψ 1 e iθ , the term in the bracket becomes Z ∂ k x ððψ 1 e iθ Þ Ã Þ∂ k y ðψ 1 e iθ Þdr − Z ∂ k y ððψ 1 e iθ Þ Ã Þ∂ k x ðψ 1 e iθ Þdr ¼ The last step is obvious because the integrated modulus of the wave function is a constant.

FIG. 1 .
FIG. 1.(a) Antidot graphene in the tight-binding model.Red and black dots represent the carbon atoms A and B. Here, L denotes the periodicity of the antidot lattice, and d denotes the diameter of the hole.(b) The matrix element magnitude plotted as a function of Bloch wave vector k.The term jH 12 j (black solid squares) is found to vary linearly with k ¼ jkj, while the intervalley scattering term jH 14 j (red solid circles) is almost a constant.(c) Band structures for various γ with a fixed periodicity L ¼ 13 nm.The black circles represent the tight-binding results, and red curves stand for the dispersion relation given by the effective Hamiltonian.The Bloch wave vector k ranges from K to Γ, then to M, and finally back to K [here, K, Γ, and M points are for the antidot (hole) reciprocal lattice].

FIG. 2 .
FIG. 2. Effective velocity v and intervalley scattering strength m are plotted as a function of the geometric factor γ. The calculated results on velocity v (in units of v F ) are shown as black solid squares, and those for m (in units of t) are shown as red open circles.These parameters are obtained for a fixed periodicity L ¼ 13 nm.The dashed curves are linear fittings.

FIG. 3 .
FIG. 3. Berry curvature matrix element Ω 11 calculated as a function of wave vector k for various values of the geometric factor γ [based on Eq. (19)].The parameter m=ℏv can be obtained from the tight-binding calculations, where the periodicity L is fixed at 13 nm.Here a denotes the carbon-carbon atomic separation (¼ 1.42 Å).

FIG. 4 .
FIG. 4. (a) Schematic illustration of an antidot graphene lattice (Lattice B) for calculating the conductance as a function of Fermi energy.The system contains 20 holes in the x direction and 10 holes in the y direction.The inset is a zoom-in image of the antidot lattice with a hole.(b) Calculated conductance plotted as a function of Fermi energy.The black curve represents the conductance of LatticeA, where the zero-conductance plateau indicates the band gap of this system.For Lattice B, the red curve exhibits no zeroconductance regions, implying that there is no band gap.Lattice C differs from Lattice B only at the edges of the holes, where we have added atomic-scale randomness.Ten different randomness profiles were calculated, and the blue curve represents the ensembleaveraged conductance.It is seen that a band gap reappears, and the magnitude of this band gap is the same as that of Lattice A, implying that the transport gap in Lattice C has the same origin as that in Lattice A, i.e., from K to K 0 intervalley scatterings.

FIG. 6 .
FIG. 6. Cartoon image illustrating different wave states' evolution under an external electric field under the influence of the Berry curvature.The value of the Berry curvature is shown in the color scale, plotted in the first Brillouin zone of the antidot graphene.The bottom left (right) shows the value of the diagonal matrix term Ω 11 (Ω 22 ), with the peak value at the Γ point of the antidot reciprocal lattice.The red wave represents the case of jη 0 1 j 2 > 1=2; it travels to the left as v x < 0. The blue wave represents the case of jη 0 1 j 2 < 1=2; it travels to the right.

FIG. 7 .
FIG. 7. Measurement setup for (a) nonlocal measurements and (b) local measurements.(c) Antidot graphene resistivity is plotted as a function of gate voltage for various temperatures (from 10 to 300 K).Inset: SEM image for the antidot graphene sample.(d) Roomtemperature conductivity as a function of gate voltage is shown to be well fitted by the Boltzmann transport theory, which gives the mobility around 1500 cm 2 =ðVsÞ and the residue carrier density n 0 ¼ 2 × 10 11 =cm 2 .

FIG. 9 .
FIG. 9. (a) Parameters f, λ and (b) coefficient αðξÞ are plotted as a function of ξ=W.The open symbols are simulation results, while the solid curves are from the analytical expression given in Ref.[11].The simulation results and the predictions of the analytical expression agree well in the range ξ ≫ W, but they differ when ξ ≤ W. The latter applies to our samples.

FIG. 11 .
FIG. 11.The measured nonlocal resistance R exp nl is plotted as a function of gate voltage for different temperatures, as shown by the solid curves.The calculated stray current R 0 nl is shown by dashed curves.The magnitude of R 0 nl is amplified by a factor of 7 for 10 K and by a factor of 4 for 20 K.

FIG. 12 .
FIG. 12.The peak values (at CNP) of measured nonlocal resistance are plotted as a function of temperature, shown as black solid squares.The calculated stray current effect R 0 nl ∝ ρðTÞ is shown by the solid red curve.The blue curve denotes the peak values of topological nonlocal resistance calculated, R nl ∝ ρ 3 ðTÞ, based on Eq. (28) for a fixed ξ ¼ 450 nm.The black dashed curve represents R nl þ R 0 nl .It is clear that there is a crossover from the high-temperature, linear ρðTÞ behavior to the low-temperature, ρ 3 ðTÞ behavior.The measured data are seen to agree with this crossover behavior (dashed curve) extremely well.

Fig. 14 .
FIG. 13.The topological nonlocal resistance (defined as R exp nl − R 0 nl ) is plotted as a function of carrier density, shown by open circles at temperatures of 10, 20, and 40 K.The fitted nonlocal resistance is shown by the solid black curves.This fitting yields a value for the decay length ξ of 450 nm.

FIG. 14 .
FIG. 14.The topological nonlocal resistance (defined as R exp nl − R 0 nl ) is plotted as a function of carrier density for sample B (upper panel), with a designed hole diameter of 80 nm, and sample C (bottom panel), with a designed hole diameter of 60 nm, shown by open circles for different samples with different geometric factors.The periodicity is fixed at 150 nm.The solid curves represent the fitted nonlocal resistance.
calculate the parameter m=ℏv (based on tight-binding calculations) and the residue carrier density n 0 (based on experimental resistivity-back gate curves).For sample B, we have m=ℏv ¼ 0.021 nm −1 and n 0 ¼ 2.2 × 10 11 =cm −2 ; for sample C, we have m=ℏv ¼ 0.015 nm −1 and n 0 ¼ 1.3 × 10 11 =cm −2 .With these input values, we fit the experimental nonlocal resistance data, shown by the black solid curves in Fig.14.These fittings yield the decay length ξ ≈ 500 nm for sample B and ξ ≈ 550 nm for sample C. It is also noted that as the diameters of the holes decrease, the signal of nonlocal resistance is reduced, due to the decrease of ρ.
In this work, we have obtained the effective Hamiltonian of antidot graphene based on the tight-binding eigenvalues and eigenfunctions.The effective Hamiltonian predicts the gap opening at the Dirac point due to intervalley scatterings.Based on the effective Hamiltonian, we find the Berry curvature to be nonzero, without breaking either the timereversal or inversion symmetry.The nonzero Berry curvature can be experimentally verified by nonlocal transport measurements, based on three separately fabricated antidot graphene samples.The very good one-parameter fitting of the nonlocal resistance data gives the topological current decay length of ξ ¼ 450, 500, and 550 nm for three samples with fixed periodicity of 150 nm and hole diameters of 100, 80, and 60 nm, respectively.ACKNOWLEDGMENTS Research support by Hong Kong Research Grant Council GRF Grants No. 16304314 and No. 16307114 is hereby gratefully acknowledged.P. S. also wishes to thank Steven Louie and C. T. Chan for helpful discussions.APPENDIX A: PROOF OF TIME-REVERSAL SYMMETRY RELATION BETWEEN DOUBLE DEGENERATE STATES By definition, under time-reversal transformation, the momentum k flips signs, and the wave function takes its complex conjugate, that is, Tψ 2 ðkÞ ¼ ψ 2 ð−kÞ Ã :