Geometric and nongeometric contributions to the surface anomalous Hall conductivity
GGeometric and nongeometric contributions to the surface anomalous Hall conductivity
Tom´aˇs Rauch, Thomas Olsen, David Vanderbilt, and Ivo Souza
1, 4 Centro de F´ısica de Materiales, Universidad del Pa´ıs Vasco (UPV/EHU), 20018 San Sebasti´an, Spain CAMD, Department of Physics, Technical University of Denmark, 2820 Kgs. Lyngby Denmark Department of Physics and Astronomy, Rutgers University, Piscataway, New Jersey 08854-8019, USA Ikerbasque Foundation, 48013 Bilbao, Spain (Dated: September 7, 2018)A static electric field generates circulating currents at the surfaces of a magnetoelectric insulator.The anomalous Hall part of the surface conductivity tensor describing such bound currents canchange by multiples of e /h depending on the insulating surface preparation, and a bulk calculationdoes not fix its quantized part. To resolve this ambiguity, we develop a formalism for calculating thefull surface anomalous Hall conductivity in a slab geometry. We identify a Berry-curvature term,closely related to the expression for the bulk anomalous Hall conductivity, whose value can changeby quantized amounts by adjusting the surface Hamiltonian. In addition, the surface anomalousHall conductivity contains a nongeometric part that does not depend on the surface preparation. I. INTRODUCTION
Certain surface properties of crystals are strongly con-strained by the bulk, and as a result they are very robustwith respect to local perturbations. An example is theareal charge density σ surf bound to an insulating surfaceof a polar insulator. For an unreconstructed defect-freesurface with outward normal ˆ n it is given by [1] σ surf = (cid:18) P + e R V c (cid:19) · ˆ n , (1)where P is the bulk electric polarization, R is a latticevector, and V c is the volume of a unit cell. Accordingto the Berry-phase theory [2], P is only defined modulo e R /V c , since it is possible to change its value by thatamount by adjusting the phases of the Bloch wave func-tions. Equation (1) assumes that a definite choice ofgauge has been made so that a unique value of P hasbeen established. (Here, the word “gauge” refers to thefreedom to adjust the phases of the Bloch eigenstates or,more generally, to perform a unitary transformation ateach k among the occupied Bloch states [3].) The secondterm in Eq. (1) amounts to an integer number of electronsper surface unit cell. Its presence is required because it isin principle possible to prepare the insulating surface indifferent ways such that the macroscopic charge per sur-face cell changes by a multiple of the elementary charge e .Thus, the quantized part of σ surf depends on the detailsat the surface but the nonquantized part does not.In this work, we consider a similar situation that arisesin insulating crystals that display the linear magnetoelec-tric (ME) effect, whereby an applied magnetic field B induces an electric polarization P , and conversely an ap-plied electric field E induces a magnetization M [4, 5].The linear ME tensor is defined as α ab = ∂P a ∂B b (cid:12)(cid:12)(cid:12)(cid:12) E =0 = ∂M b ∂ E a (cid:12)(cid:12)(cid:12)(cid:12) B =0 . (2)The full ME response contains both frozen-ion andlattice-mediated contributions, and each can be further decomposed into spin and orbital parts. In the following,we focus exclusively on the frozen-ion orbital response.The bulk magnetization generated by a static electricfield gives rise to surface currents K = M × ˆ n . (3)In the case of an insulating surface, this is the full cur-rent response. It is described at linear order by a 2 × σ surf ab = ∂K a /∂ E b , and thesurface anomalous Hall conductivity (AHC) is defined asthe antisymmetric part of the 2 × σ AHsurf ˆ n where σ AHsurf = − (cid:15) cdb σ surf cd ˆ n b , (4)the surface anomalous Hall current density becomes K AH = σ AHsurf ˆ n × E . (5)From Eqs. (2) and (3) we find σ surf cd = ∂K c ∂ E d = ∂∂ E d (cid:15) cea M e ˆ n a = (cid:15) cea α de ˆ n a , (6)and plugging this expression into Eq. (4) leads to σ AHsurf = −
12 Tr( α ) + 12 α ab ˆ n a ˆ n b . (7)Separating the ME tensor on the right-hand side into anisotropic trace piece and a traceless part, α ab = α iso δ ab + (cid:101) α ab , (8)we arrive at the relation σ AHsurf := − α iso + 12 (cid:101) α ab ˆ n a ˆ n b . (9)We use the special symbol := to indicate that while theleft-hand side is uniquely defined for a given surface ter-mination, the right-hand side carries a quantum of inde-terminacy, since the bulk quantity α iso is gauge invariantonly modulo e /h [6, 7]. a r X i v : . [ c ond - m a t . m e s - h a ll ] S e p Once a definite value has been chosen for the multival-ued quantity α iso , Eq. (9) can be rewritten as σ AHsurf = − α iso + m e h + 12 (cid:101) α ab ˆ n a ˆ n b . (10)Different choices of the integer m in the middle termcorrespond to different surface preparations exhibitingvalues of the surface AHC that differ by a multiple of thequantum of conductance. This integer can be changed,in principle, by stitching a quantum anomalous Hall layerto the surface [7, 8] or otherwise changing the surfaceHamiltonian, or by means of an adiabatic pumping cyclecharacterized by a nonzero second Chern number [10, 11].It follows that only the nonquantized part of the surfaceAHC is a bulk property, in close analogy with Eq. (1)for the surface charge. These features are described bythe phenomenology of axion electrodynamics [6, 12], and α iso is sometimes referred to as the axion ME coupling.We are now ready to formulate the main question be-hind this work. Suppose we have a ME insulator (itshould break both inversion and time reversal symme-try), and we consider a specific insulating surface. Howcan we calculate the surface AHC, not just up to a quan-tum but exactly? Since we are given a definite surfaceHamiltonian, there should be a definite answer withoutany quantum of ambiguity. We shall answer this questionby developing a formalism that allows one to calculate thesurface AHC unambiguoulsy using a slab geometry.The paper is organized as follows. We begin in Sec. IIby calculating, at linear order, the local current responseof an insulator to a static homogeneous electric field. Thelocal AHC, defined as the antisymmetric part of this lo-cal conductivity tensor, is then separated into geometricand nongeometric parts. Starting from the expression forthe local AHC, we obtain in Sec. III an expression for thesurface AHC of a slab, which we again separate into ge-ometric and nongeometric parts. In Sec. IV we calculatenumerically the surface AHC for slabs of tight-binding(TB) models and compare the results, via Eq. (10), withindependent calculations of the bulk ME tensor. We con-clude in Sec. V with a summary, and leave a lengthierderivation to an appendix. There are two scenarios compatible with Eq. (10). If the integer m is the same for all crystal facets, the entire surface is insulatingand the term me /h gives an isotropic contribution to the surfaceAHC [8]. If adjacent facets have different m values, there arechiral conducting channels along the connecting hinges [9]. II. LOCAL ANOMALOUS HALLCONDUCTIVITYA. Linear-response calculation
The local conductivity and local AHC are defined as σ ab ( r ) = ∂j a ( r ) ∂ E b (cid:12)(cid:12)(cid:12)(cid:12) E =0 (11)and σ AH c ( r ) = − (cid:15) abc σ ab ( r ) , (12)respectively. E denotes a static homogeneous electricfield and j ( r ) is the microscopic induced current den-sity, whose anomalous Hall part is given by j AH ( r ) = σ AH ( r ) × E .We wish to calculate the local AHC for an insulat-ing medium at zero temperature described by a single-particle Hamiltonian ˆ H . The current operator isˆ j ( r ) = − e | r (cid:105)(cid:104) r | ˆ v + ˆ v | r (cid:105)(cid:104) r | ) , (13)where ˆ v = (1 /i (cid:126) )[ˆ r , ˆ H ] is the velocity operator and e > j ( r ) = Tr (cid:104) ˆ P ˆ j ( r ) (cid:105) = − e Re (cid:104) r | ˆ v ˆ P | r (cid:105) , (14)where ˆ P denotes the projection operator onto the occu-pied states. An expression for the local conductivity (11)can now be obtained by differentiating Eq. (14) with re-spect to E . For that purpose we write ˆ H = ˆ H + e E · ˆ r where ˆ H is the unperturbed Hamiltonian, and note thatsince [ˆ r a , ˆ r b ] = 0 the operator ˆ v reduces to (1 /i (cid:126) )[ r , ˆ H ].Hence the electric field enters Eq. (14) via ˆ P only, lead-ing to σ ab ( r ) = ( − e )Re (cid:104) r | ˆ v a ∂ E b ˆ P | r (cid:105) , and inserting thisexpression in Eq. (12) we arrive at σ AH ( r ) = e (cid:104) r | ˆ v × ∂ E ˆ P | r (cid:105) . (15)Finally, from first-order perturbation theory we get ∂ E ˆ P = − e (cid:88) v,c (cid:18) | c (cid:105) (cid:104) c | ˆ r | v (cid:105) E cv (cid:104) v | + | v (cid:105) (cid:104) v | ˆ r | c (cid:105) E cv (cid:104) c | (cid:19) , (16)where | v (cid:105) and | c (cid:105) denote occupied and empty energyeigenstates, respectively, and E cv = E c − E v . Equa-tions (15) and (16) give the full local AHC; below, weseparate it into geometric and nongeometric parts. We are ignoring local-field corrections, which introduce a depen-dence of ˆ H on E through the self-consistent charge density. Suchterms are not difficult to derive, but they are absent from ournon-self-consistent TB calculations. B. Separation of the local AHC into geometric andnongeometric parts
Consider the isotropic ME response of a bounded sam-ple, defined as a iso = 13 (cid:88) a =1 ∂m a ∂ E a (17)in terms of the orbital moment m = 12 (cid:90) r × j ( r ) d r. (18)For a globally insulating crystallite of volume V , a iso /V converges in the V → ∞ limit to one of the multiplevalues of α iso , with the specific value depending on thesurface preparation [8]. Plugging Eq. (18) into Eq. (17)and comparing with the definition of the local AHC inEq. (12), we find a iso = − (cid:90) r · σ AH ( r ) d r. (19)This relation will be used below to isolate the geometricand nongeometric contributions to the local AHC, butfirst we need some results from the microscopic theory ofthe orbital ME response in insulators [13, 14].The quantum-mechanical expression for the bulk MEtensor α ab comprises an isotropic geometric term knownas the Chern-Simons (CS) term, and a nongeometric termknown as the Kubo or cross-gap (cg) term that has bothisotropic and anisotropic parts [13, 14]. The relation be-tween those two terms and the decomposition in Eq. (8)can be summarized as follows: α ab = α cg ab (cid:122) (cid:125)(cid:124) (cid:123)(cid:0) α CS + α cgiso (cid:1) δ ab + (cid:101) α ab (cid:124) (cid:123)(cid:122) (cid:125) α iso . (20)The expressions for α CS and α cg ab take the form of integralsover the Brillouin zone (BZ), and can be found in Refs. 13and 14. In the case of α CS , the integrand only containsthe unperturbed cell-periodic Bloch functions and theirfirst k derivatives; it is a gauge-dependent quantity, butafter integration over the entire BZ it becomes gauge in-variant modulo e /h . The expression for α cg ab contains inaddition the perturbed wave functions and the velocityoperator, and is fully gauge invariant at each point in theBZ.The ME tensor a ab = ( ∂m b /∂ E a ) B =0 of a finite crys-tallite can be similarly decomposed into geometric andnongeometric terms [13]. Because surface contributionsare now included, the geometric part a CS of the isotropicpiece a iso is unique for a given surface preparation. It isgiven by [13] a CS = − πe h (cid:15) abc Im Tr (cid:104) ˆ P ˆ r a ˆ P ˆ r b ˆ P ˆ r c (cid:105) = 2 πe h (cid:15) abc (cid:90) r c Im (cid:104) r | ˆ P ˆ r a ˆ Q ˆ r b ˆ P | r (cid:105) d r, (21) where ˆ P = (cid:80) v | v (cid:105)(cid:104) v | and ˆ Q = ˆ − ˆ P are the ground-state projector and its complement, respectively. Com-paring with Eq. (19), we are led to identify a geometric(CS) contribution to the local AHC given by σ AHCS ( r ) = e h C ( r ) , (22a) C ( r ) = − π Im (cid:104) r | ˆ P ˆ r ˆ Q × ˆ Q ˆ r ˆ P | r (cid:105) . (22b)One can obtain the nongeometric (cross-gap) part ofthe local AHC in a similar way, starting from the nonge-ometric part of the orbital ME coupling of a crystallite.This is done in the Appendix and as expected the resultis that the cross-gap local AHC is equal to the differencebetween the full local AHC (15) and the CS term above, σ AHcg ( r ) = σ AH ( r ) − σ AHCS ( r ) . (23)Equations (15), (22), and (23) are the main results ofthis section. C. Discussion
The appearance of a nongeometric term in the localAHC may seem surprising at first, since the intrinsicmacroscopic AHC of a bulk crystal or slab is known tobe purely geometric: it is given by the BZ integral of theBerry curvature of the occupied Bloch states [15]. Theexplanation, as we will see Sec. III D, is that the non-geometric term always integrates to zero across the fullwidth of a slab, dropping out from the net AHC of theslab. As will become clear in the following, that termdoes contribute to the AHC of a single surface.The nongeometric part of the local AHC was over-looked in some previous studies [7, 16], where the localAHC was formulated as a spatially-resolved Berry curva-ture. As for the geometric part, the expression in Eq. (22)is consistent with the previous literature. Consider a flatcrystallite lying on the ( x, y ) plane, and integrate thequantity C z ( r ) given by Eq. (22b) over all z to obtaina dimensionless quantity C z ( x, y ). This is precisely the“local Chern number” introduced in Ref. [17]. For a slab,the average of Eq. (22a) over a two-dimensional (2D) cellat fixed z is essentially identical to the “layer-resolvedAHC” of Ref. [7].In the next section we calculate the layer-resolvedAHC including both geometric and nongeometric con-tributions, and use it to evaluate the surface AHC. III. SURFACE ANOMALOUS HALLCONDUCTIVITYA. Evaluation in a slab geometry
Consider an insulating slab with the outward normalˆ n = ˆ z of the top surface pointing along a reciprocal-lattice vector b . We assume that the slab thickness L ismuch larger than the lattice constant c = 2 π/ | b | in thesurface-normal direction. We also assume a defect-freesurface, and introduce a “layer-resolved” AHC for theslab by averaging the z component of the local AHC (12)over a surface unit cell at fixed z : σ AHslab ( z ) = 1 A c (cid:90) A c σ AH z ( x, y, z ) dx dy. (24)The net AHC of a slab is given by σ AHslab = (cid:90) σ AHslab ( z ) dz, (25)where the range of integration is chosen to span the fullwidth of the slab, including the exponential tails of thewave functions outside the two surface regions. For an in-sulating slab the result is quantized in units of e /h [15]: σ AHslab = e h C slab , (26)where the integer C slab is the Chern number of the slab(see Sec. III D). In order for this equation to be mean-ingful, we are assuming that the slab is cut from a bulkinsulator for which all three of the bulk Chern indices C j [18] are zero.As a first step towards calculating the surface AHC,we filter out the atomic-scale oscillations in the layer-resolved AHC by performing a “sliding-window average”over one vertical lattice constant: σ AHslab ( z ) = 1 c (cid:90) z + c/ z − c/ σ AHslab ( z (cid:48) ) dz (cid:48) . (27)Because we assumed C = 0, this coarse-grained AHCmust vanish exponentially in the bulklike interior regionof the slab, and it can only become nonzero near the twosurfaces. The macroscopic AHC of the top surface cannow be expressed as σ AHsurf = (cid:90) z σ AHslab ( z ) dz, (28)with z chosen in the bulklike region of the slab, and theupper limit of integration placed at an arbitrary point in the vacuum region above the top surface. The AHC ofthe bottom surface is ( e /h ) C slab − σ AHsurf .For numerical work, it is more convenient to recastEq. (28) as σ AHsurf = (cid:90) σ AHslab ( z ) f ramp ( z − z ) dz, (29)where the range of integration spans the full width of theslab, and f ramp is a ramp-up function defined as f ramp ( z ) = , for z < − c/ z/c + 1 / , for − c/ < z < c/ , for z > c/ . (30)To summarize, Eq. (29) gives the surface AHC interms of the layer-resolved AHC of Eq. (24), for whichwe provide an explicit formula below. B. Layer-resolved anomalous Hall conductivity
We evaluate the layer-resolved AHC by insertingEq. (15) for the local AHC into Eq. (24). The ground-state projector expressed in terms of the valence eigen-states of the slab reads asˆ P = 1 N (cid:88) k v | ψ k v (cid:105)(cid:104) ψ k v | (31)(the summation in k = ( k x , k y ) is over a uniform mesh of N points covering the surface BZ), and we need its linearchange under an in-plane electric field, ∂ E ˆ P = 1 N (cid:88) k v e i k · ˆ r (cid:16) | (cid:101) ∂ E u k v (cid:105)(cid:104) u k v | + | u k v (cid:105)(cid:104) (cid:101) ∂u k v | (cid:17) e − i k · ˆ r . (32)Here, | u k v (cid:105) is the cell-periodic part of | ψ k v (cid:105) , and | (cid:101) ∂ E u k v (cid:105) = ie (cid:88) c | u k c (cid:105) (cid:126) v k cv E k cv (33)is the projection of | ∂ E u k v (cid:105) onto the conduction bands,with v k cv = (cid:104) u k c | ˆ v k | u k v (cid:105) and ˆ v k = e − i k · ˆ r ˆ v e i k · ˆ r . Equa-tion (32) is essentially the same as Eq. (16), but adaptedto a slab geometry. Putting everything together and let-ting N → ∞ , we arrive at σ AHslab ( z ) = e πh (cid:90) d k (cid:90) A c dx dy (cid:88) v Re (cid:104) (cid:104) r | (cid:126) ˆ v k × (cid:16) | (cid:101) ∂ E u k v (cid:105)(cid:104) u k v | r (cid:105) + | u k v (cid:105)(cid:104) (cid:101) ∂ E u k v | r (cid:105) (cid:17)(cid:105) z , (34)where the integral in k is over the surface BZ. C. Separation of the layer-resolved AHC intogeometric and nongeometric parts
To find the geometric part of the layer-resolved AHC,we repeat the above steps, simply replacing Eq. (15) forthe full local AHC with Eq. (22) for the geometric part.Inserting Eq. (31) for ˆ P in Eq. (22b) gives C z ( r ) = − πN Im (cid:88) kk (cid:48) vv (cid:48) ψ ∗ k v ( r ) ψ k (cid:48) v (cid:48) ( r ) (cid:104) ψ k v (cid:48) | ˆ x ˆ Q ˆ y | ψ k v (cid:105) . (35)Writing ˆ Q as (1 /N ) (cid:80) k c | ψ k c (cid:105)(cid:104) ψ k c | and using the iden-tity (cid:104) ψ k v | ˆ r j | ψ k (cid:48) c (cid:105) = iN (cid:104) u k v | ∂ k j u k c (cid:105) δ kk (cid:48) valid for off-diagonal position matrix elements [19], we obtain C z ( r ) = − πN Im (cid:88) k u ∗ k v ( r ) u k v (cid:48) ( r ) F xy k v (cid:48) v , (36)where F xy k v (cid:48) v = (cid:88) c (cid:104) ∂ k x u k v (cid:48) | u k c (cid:105)(cid:104) u k c | ∂ k y u k v (cid:105) = ( F yx k vv (cid:48) ) ∗ (37)is the metric-curvature tensor [3]. Inserting Eq. (36)in Eq. (22a) for σ AHCS ( r ) and plugging the result intoEq. (24) for the layer-resolved AHC yields, for N → ∞ , σ AHslab , CS ( z ) = − e πh (cid:90) Im Tr [ O k ( z ) F xy k ] d k. (38)The trace is over the valence bands, and O k vv (cid:48) ( z ) = (cid:90) A c u ∗ k v ( x, y, z ) u k v (cid:48) ( x, y, z ) dx dy (39)is a layer-resolved overlap matrix.Equation (38) can be brought to a more transparentform by decomposing the metric-curvature tensor intoHermitian and anti-Hermitian parts in the band indicesas F xy k nm = R xy k nm + (1 / i ) (cid:101) Ω xy k nm , where R xy k nm = 12 F xy k nm + 12 ( F xy k mn ) ∗ (40)is the quantum metric and (cid:101) Ω xy k nm = iF xy k nm + ( iF xy k mn ) ∗ (41)is the gauge-covariant Berry curvature, related to theBerry connection A a k nm = i (cid:104) u k n | ∂ k a u k m (cid:105) and to the non-covariant Berry curvature Ω xy k nm = ∂ k x A y k nm − ∂ k y A x k nm by (cid:101) Ω xy k nm = Ω xy k nm − i [ A x k , A y k ] nm . (42)Since the matrices O k ( z ), R xy k , and (cid:101) Ω xy k are all Hermitianand the trace of the product of two Hermitian matricesis real, R xy k drops out from Eq. (38), which reduces to σ AHslab , CS ( z ) = e πh (cid:90) Tr (cid:104) O k ( z ) (cid:101) Ω xy k (cid:105) d k. (43)This expression for the CS layer-resolved AHC in termsof the layer-resolved overlap matrix and the non-AbelianBerry curvature is the central result of this section. Equation (43), which essentially agrees with Eq. (12)of Ref. [7], only accounts for part of the layer-resolvedAHC, whose full amount is given by Eq. (34). Accordingto Eq. (23) for the local AHC, the remainder is the non-geometric (or cross-gap) part of the layer-resolved AHC, σ AHslab , cg ( z ) = σ AHslab ( z ) − σ AHslab , CS ( z ) . (44)To review, the surface AHC is calculated from Eq. (29),where we insert either Eq. (34) to obtain the grand total,or Eq. (43) to find the geometric part. D. Discussion
It was already mentioned in Sec. II C that although thenongeometric part of the layer-resolved AHC can give anet contribution to the surface AHC, its contribution tothe AHC of the entire slab vanishes. To establish thisresult, let us show that the geometric part of the slabAHC coincides with the total.We begin with the total slab AHC. Inserting Eq. (34)for the layer-resolved AHC into Eq. (25) yields σ AHslab = − e πh (cid:90) d k Im (cid:88) vc (cid:126) v x k vc v y k cv E k cv − ( x ↔ y ) . (45)Using (cid:126) v k vc = − iE k cv A k vc to remove the energy denom-inator and noting that Im (cid:80) k vv (cid:48) A x k vv (cid:48) A y k v (cid:48) v = 0, thesum over conduction bands c can be replaced by a sumover all bands n (the term n = v vanishes). Comparingwith the Berry curvature Ω xy k v = − (cid:80) n (cid:54) = v A x k vn A y k nv we arrive at Eq. (26) for the total slab AHC, with theslab Chern number given by [15] C slab = 12 π (cid:90) d k (cid:88) v Ω xy k v . (46)To obtain the gometric part of the slab AHC, insertEq. (42) into Eq. (43) and plug the result into Eq. (25).Using (cid:82) O k vv (cid:48) ( z ) dz = δ vv (cid:48) we get σ AHslab , CS = C slab ( e /h ),which is the same as the total slab AHC. Thus, Eq. (44)must integrate to zero across the entire slab.The net amount of AHC contributed by the cross-gapterm (44) across a single surface region is equal to σ AHsurf , cg = − α cgiso + 12 (cid:101) α ab ˆ n a ˆ n b , (47)and the CS term (43) contributes the additional amount σ AHsurf , CS = − α CS + m e h . (48)Together, they make up the full surface AHC of Eq. (10). To obtain the expression in Ref. [7] for the CS layer-resolvedAHC starting from the CS local AHC, one can repeat the deriva-tion of Eq. (43) with a single modification: in Eq. (22b) for C ( r ),exchange ˆ P and ˆ Q and remove the minus sign. It can be eas-ily verified that this “particle-hole transformation” leaves C ( r )unchanged, as expected on physical grounds. IV. NUMERICAL RESULTS
In the following, we present the results of slab calcu-lations of the surface AHC for two different TB models,and compare the results with bulk calculations of the or-bital ME tensor [13]. First, let us briefly describe our TBimplementation of the surface AHC expressions.
A. Tight-binding formulation of the surface AHC
In the TB context, the integration over z in Eq. (29)for the surface AHC gets replaced by a summation overa layer index l . The ramp-up function is evaluated atthe discrete layer coordinates z ( l ), and the layer-resolvedAHC becomes σ AHslab ( l ).In Eq. (34) for the layer-resolved AHC, | r (cid:105) is replacedby | i (cid:105) representing a TB orbital φ i ( r ) = (cid:104) r | i (cid:105) , the inte-gration (cid:82) A c dx dy is replaced by a summation (cid:80) i ∈ l overthe orbitals within one surface unit cell in layer l , and thematrix elements of the velocity operator are evaluated bymaking the diagonal approximation (cid:104) i | ˆ r | j (cid:105) = τ i δ ij for theposition operator in the TB basis [20].In Eq. (43) for the CS layer-resolved AHC the overlapmatrix becomes O k vv (cid:48) ( l ) = (cid:88) i ∈ l u ∗ k v ( i ) u k v (cid:48) ( i ) , (49)and the non-Abelian Berry curvature can be evaluatedfrom Eqs. (37) and (41) with the help of the identity (cid:104) u k c | ∂ k u k v (cid:105) = − (cid:104) u k c | (cid:126) ˆ v k | u k v (cid:105) E k cv . (50) B. Anisotropic cubic-lattice model
As our first test case, we consider a model of a ME insu-lator with no symmetry. This provides the most challeng-ing case for the theory, since all nine components of theME tensor are nonzero and different from one another.The resulting surface AHC has both CS and cross-gapcontributions, and it varies with the surface orientation.We choose the TB model described in Appendix Aof Ref. [13]. This is a spinless model defined on a cu-bic lattice, with one orbital per site and eight sites percell, where inversion and time-reversal symmetry are bro-ken by assigning random on-site energies and complexfirst-neighbor hoppings of fixed magnitude. We take themodel parameters listed in Table A.1 of Ref. [13], andchoose the two lowest bands to be the valence bands. Asin that work, all parameters are kept fixed except for onehopping phase ϕ which is scanned from 0 to 2 π , and theresults are plotted as a function of this phase ϕ .This model is intended as a model for a conventionalME insulator, in which the isotropic response a iso /V of aninsulating crystallite is very small relative to the quantum layer l -0.4-0.20.00.20.4 σ AH s l a b ( × − e h ) totalCS FIG. 1. Coarse-grained layer-resolved AHC [Eq. (51)] for a16-atom-thick slab of the cubic-lattice model with ϕ = π . e /h . The surface AHC of such a system is most natu-rally described by setting m = 0 in Eqs. (10) and (48)while choosing α iso and α CS in the range [ − e / h, e / h ].(In the next subsection, we will consider a model withthe opposite characteristics, i.e., with a large isotropicME response of the order of e /h .)We construct a slab with a thickness of 16 atomic layers(8 lattice constants) along z . The layer-resolved AHCdisplays strong oscillations from one layer to the next,which we filter out by averaging over two consecutivelayers: σ AHslab ( l + 1 /
2) = 12 (cid:2) σ AHslab ( l ) + σ AHslab ( l + 1) (cid:3) . (51)This quantity is plotted as the solid line in Fig. 1 for ϕ = π , and the dashed line shows the CS contributionthe CS layer-resolved AHC was calculated for a differentTB model in Ref. [7], and is shown in Fig. 2 therein).Both quantities are nonzero in the surface regions only,quickly dropping to almost zero within four subsurfacelayers. The fact that the two curves in Fig. 1 are not per-fectly odd about the center of the slab can be attributedto the lack of mirror symmetry in the model. We havechecked that both (cid:80) l σ AHslab ( l ) and (cid:80) l σ AHslab , CS ( l ) vanishidentically, as should be the case for a slab with C slab = 0,so that the macroscopic surface AHC is equal and oppo-site for the two surfaces. On a given surface, the CS partof the AHC has the opposite sign compared to the to-tal. The cross-gap contribution therefore prevails, as inordinary ME insulators [8].The macroscopic AHC of the top surface is plotted ver-sus ϕ as the solid line in the top panel of Fig. 2, wherethe CS contribution is again shown as a dashed line. Forcomparison, we plot as filled (total) and empty (CS) cir-cles the quantities on the right-hand side of Eqs. (10)and (48), respectively (with m = 0). The precise agree-ment validates our expression for the surface AHC, andits decomposition into geometric and nongeometric parts. -0.9-0.6-0.30.00.3 ˆ n = ˆ z σ AH s u r f ( × − e h ) π/ π π/ π ϕ -0.9-0.6-0.30.00.3 ˆ n = ˆ x FIG. 2. AHC of the top surface (upper panel) and of theright surface (lower panel) of 16-layer slabs of the cubic-latticemodel cut along z and x , respectively, as a function of thecyclic parameter ϕ . The solid (dashed) lines denote the total(CS) surface AHC. Circles represent the quantity appearingon the right-hand side of Eq. (9), with filled and empty circlesdenoting the total and the CS piece, respectively. In the lower panel of Fig. 2 we show results for a slabcut along x . The total surface AHC is different from thaton the upper panel, as expected for an anisotropic modelfrom the last term in Eq. (47). The CS surface AHC isthe same in both panels, confirming that its nonquan-tized part does not depend on the surface orientation, asexpected from Eq. (48). C. Layered Haldane model
We now turn to a model for which α iso and α CS arenot always small compared to the quantum e /h , so thatthe branch choice in Eqs. (10) and (48) becomes ambigu-ous. We choose the TB model introduced in Ref. [11].This is a layered model on a hexagonal lattice, with fourorbitals per cell. It can be obtained by stacking Haldane-model [21] layers with alternating parameters, and thencoupling them via interlayer hopping terms. The modelacts as a quantum pump of axion ME coupling: a slow π/ π π/ π φ σ AH s u r f ( e h ) FIG. 3. Cyclic evolution of the Hamiltonian of a 10-layer slabof the layered Haldane model, during which the isotropic MEcoupling α CS of the bulk crystal changes by e /h . The AHCof the top surface is plotted as a solid line, and the blackand gray circles denote two different branches of − α CS . Fora given choice of branch, Eq. (48) is satisfied throughout thecycle with the value of m increasing by one at φ c = 3 π/ periodic variation in its parameters can gradually change α iso by e /h over one cycle. Below, we monitor the evo-lution of the surface AHC during one pumping cycle.We begin by noting that the cross-gap contribution tothe ME tensor vanishes identically for this model. Thereason can be found in Ref. [14], where a set of condi-tions were derived under which α cg ab vanishes in certainfour-band models having some kind of particle-hole sym-metry. Those conditions hold for several models proposedin the literature including the present one, so that the to-tal surface AHC (10) reduces to the CS part (48).We construct a slab containing 10 hexagonal layersstacked along z . The pumping cycle is parametrized byan angle φ that modulates the model parameters accord-ing to Eqs. (57c) and (61) in Ref. [11]. As in Sec. III.Dof that work, we consider the situation where the entireslab, including the surfaces, returns to its initial state atthe end of the cycle,ˆ H slab ( φ = 2 π ) = ˆ H slab ( φ = 0) . (52)The quantities σ AHsurf ( φ ) and − α CS ( φ ) are plotted inFig. 3 as solid lines and filled circles, respectively; the lat-ter is a multivalued quantity, and two different branchesare shown as black and gray circles. Equation (48) as-sumes that a specific branch has been chosen, and wepick the one represented by the black circles. With thatchoice, Eq. (48) is satisfied with m = 0 for 0 < φ < π/ m = 1 for 3 π/ < φ < π .The actual value of m at each φ depends on the par-ticular gauge choice, but the important point is thatEq. (48) cannot be satisfied keeping m fixed for all φ .Equation (52) implies σ AHsurf ( φ = 2 π ) = σ AHsurf ( φ = 0), andthe only way this can be reconciled with the pumping of π/ π π/ π φ (surface layers) σ AH s u r f ( e h ) FIG. 4. The same quantities as in Fig. 3, but now for acyclic evolution of the Hamiltonian of the surface layers only,keeping the Hamiltonian of the rest of the slab fixed.
CS axion coupling in the bulk region, α CS ( φ = 2 π ) = α CS ( φ = 0) + e h , (53)is if the integer m in Eq. (48) increases by one during thecycle. This change in the quantized part of the surfaceAHC is caused by a topological phase transition at thesurface: the energy gap of the surface bands closes at φ c = 3 π/
2, forming a Weyl point in ( k x , k y , φ )-space thattransfers a quantum of Berry-curvature flux between thevalence and conduction bands as φ crosses φ c [11].It is also possible to change the quantized part of thesurface AHC by adjusting only the surface Hamiltonian.This is illustrated in Fig. 4, where we held the Hamilto-nian of all subsurface layers fixed at φ = 0, but modu-lated the onsite energy of the top and bottom layers by φ according to Eq. (57c) in Ref. [11]. Now the bulk MEcoupling α CS is held at zero for all φ , as indicated by theblack circles. The surface AHC vanishes during half ofthe cycle leading to m = 0 in Eq. (48), and it becomes − e /h during the other half where m = −
1. At the crit-ical points φ c = π/ φ c = 3 π/
2, the surface statesbecome gapless.
V. SUMMARY
In this work we have derived practical expressions forcalculating the full surface AHC of an insulating slab,including the quantized part that depends on the surfacepreparation. That quantized part resides in a geometricterm written in terms of the gauge-covariant Berry cur-vature matrix of the slab wave functions. The full surfaceAHC contains an additional nongeometric term. Like thenonquantized part of the geometric surface AHC, thatterm is only apparently a surface property, but is in factfully determined by the bulk ME tensor [13, 14]. Numerical TB calculations were carried out to showthat our expressions satisfy the phenomenological rela-tion in Eq. (10) between the surface AHC and the bulkME coupling. The ability to change the surface AHC bymultiples of e /h by adjusting the surface Hamiltonianwas illustrated for a layered Haldane model.The formalism developed in this work provides a sim-pler way of determining the quantized part of the sur-face AHC than an alternative approach based on hybridWannier functions [11]. It could be particularly useful forcharacterizing the nontrivial surface topology in second-order three-dimensional topological insulators with chiralhinge states [9, 22]. ACKNOWLEDGMENTS
This work was supported by the ForschungsstipendiumGrant No. RA 3025/1-1 from the Deutsche Forschungs-gemeinschaft (T. R.), by Grant No. FIS2016-77188-Pfrom the Spanish Ministerio de Econom´ıa y Competi-tividad (T. R. and I. S.), and by NSF Grant No. DMR-1629059 (D. V.).
APPENDIX: DERIVATION OF EQ. (23) FORTHE CROSS-GAP LOCAL AHC
The orbital moment (18) of a finite sample in a staticelectric field can be decomposed as [13] m ( E ) = m CS ( E ) + m cg ( E ) , (A1a) m CS ( E ) = − πe h (cid:15) ijl Im Tr (cid:104) ˆ P ˆ r i ˆ P ˆ r j ˆ P ˆ r l (cid:105) E , (A1b) m cg i ( E ) = πeh (cid:15) ijl Im Tr (cid:104) ˆ P ˆ r j ˆ Q ˆ H ˆ Q ˆ r l − ˆ Q ˆ r j ˆ P ˆ H ˆ P ˆ r l (cid:105) . (A1c)As in Sec. II A, ˆ P denotes the projection operator ontothe occupied states in the presence of the field, ˆ Q = ˆ − ˆ P ,and ˆ H is the unperturbed Hamiltonian. Note that theCS term has an explicit linear dependence on E , whilethe cross-gap (cg) term only depends on E implicitly viathe projection operators. With the above decomposition,the isotropic ME coupling (17) becomes a iso = a CSiso + a cgiso , (A2)with a CSiso given by Eq. (21) and a cgiso = − (cid:90) r i (cid:15) ijl ∂T j ( r , E ) ∂ E l (cid:12)(cid:12)(cid:12)(cid:12) E =0 d r, (A3)where T j ( r , E ) = πeh Im (cid:104) r | ˆ P ˆ r j ˆ Q ˆ H ˆ Q − ˆ Q ˆ r j ˆ P ˆ H ˆ P | r (cid:105) . (A4)Comparing Eq. (A3) with Eq. (19) for a iso we concludethat the cross-gap local AHC is given by σ AHcg ,i ( r ) = (cid:15) ijl ∂T j ( r , E ) ∂ E l (cid:12)(cid:12)(cid:12)(cid:12) E =0 . (A5)Our remaining task is to show that this expression isequivalent to Eq. (23). We begin by plugging Eq. (A4) into Eq. (A5). Thisgenerates a total of six terms, σ AHcg ,i ( r ) = πeh (cid:15) ijl Im (cid:104) r | (cid:104) − ˆ P ˆ r j (cid:16) ˆ H ˆ Q ˆ P l (cid:17) − ˆ P ˆ r j (cid:16) ˆ P l ˆ Q ˆ H (cid:17) + ˆ P l ˆ r j ˆ H ˆ Q − ˆ Q ˆ r j (cid:16) ˆ H ˆ P ˆ P l (cid:17) − ˆ Q ˆ r j (cid:16) ˆ P l ˆ P ˆ H (cid:17) + ˆ P l ˆ r j ˆ H ˆ P (cid:105) | r (cid:105) , (A6)where we used the fact that ˆ H commutes with both ˆ P and ˆ Q , and introduced the notation ˆ P l = ∂ E l ˆ P = − ∂ E l ˆ Q for the Cartesian components of ∂ E ˆ P . Using Eq. (16) forthat operator, the individual terms in Eq. (A6) become − Im (cid:104) r | ˆ P ˆ r j (cid:16) ˆ H ˆ Q ˆ P l (cid:17) | r (cid:105) = eE c E cv Im [ (cid:104) r | v (cid:48) (cid:105)(cid:104) v (cid:48) | ˆ r j | c (cid:105)(cid:104) c | ˆ r l | v (cid:105)(cid:104) v | r (cid:105) ] , (A7a) − Im (cid:104) r | ˆ P ˆ r j (cid:16) ˆ P l ˆ Q ˆ H (cid:17) | r (cid:105) = eE c E cv Im [ (cid:104) r | v (cid:48) (cid:105)(cid:104) v (cid:48) | ˆ r j | v (cid:105)(cid:104) v | ˆ r l | c (cid:105)(cid:104) c | r (cid:105) ] , (A7b)Im (cid:104) r | ˆ P l ˆ r j ˆ H ˆ Q | r (cid:105) = − eE c (cid:48) E cv Im [ (cid:104) r | c (cid:105)(cid:104) c | ˆ r l | v (cid:105)(cid:104) v | ˆ r j | c (cid:48) (cid:105)(cid:104) c (cid:48) | r (cid:105) ] − eE c (cid:48) E cv Im [ (cid:104) r | v (cid:105)(cid:104) v | ˆ r l | c (cid:105)(cid:104) c | ˆ r j | c (cid:48) (cid:105)(cid:104) c (cid:48) | r (cid:105) ] , (A7c) − Im (cid:104) r | ˆ Q ˆ r j (cid:16) ˆ H ˆ P ˆ P l (cid:17) | r (cid:105) = eE v E cv Im [ (cid:104) r | c (cid:48) (cid:105)(cid:104) c (cid:48) | ˆ r j | v (cid:105)(cid:104) v | ˆ r l | c (cid:105)(cid:104) c | r (cid:105) ] , (A7d) − Im (cid:104) r | ˆ Q ˆ r j (cid:16) ˆ P l ˆ P ˆ H (cid:17) | r (cid:105) = eE v E cv Im [ (cid:104) r | c (cid:48) (cid:105)(cid:104) c (cid:48) | ˆ r j | c (cid:105)(cid:104) c | ˆ r l | v (cid:105)(cid:104) v | r (cid:105) ] , (A7e)Im (cid:104) r | ˆ P l ˆ r j ˆ H ˆ P | r (cid:105) = − eE v (cid:48) E cv Im [ (cid:104) r | c (cid:105)(cid:104) c | ˆ r l | v (cid:105)(cid:104) v | ˆ r j | v (cid:48) (cid:105)(cid:104) v (cid:48) | r (cid:105) ] − eE v (cid:48) E cv Im [ (cid:104) r | v (cid:105)(cid:104) v | ˆ r l | c (cid:105)(cid:104) c | ˆ r j | v (cid:48) (cid:105)(cid:104) v (cid:48) | r (cid:105) ] , (A7f)where a summation over repeated band indices is im-plied. We wish to bring the sum of all these terms into a“cross-gap” form, where dipole matrix elements only con-nect occupied and empty states. Four of the eight termsabove already have that form and they can be combinedin pairs, (A7a) with the second term in (A7f) and thefirst term in (A7c) with (A7d), to get e ( E c + E v (cid:48) ) E cv Im [ (cid:104) r | v (cid:48) (cid:105)(cid:104) v (cid:48) | ˆ r j | c (cid:105)(cid:104) c | ˆ r l | v (cid:105)(cid:104) v | r (cid:105) ]+ e ( E v + E c (cid:48) ) E cv Im [ (cid:104) r | c (cid:48) (cid:105)(cid:104) c (cid:48) | ˆ r j | v (cid:105)(cid:104) v | ˆ r l | c (cid:105)(cid:104) c | r (cid:105) ] . (A8)In the remaining four terms, we use the completenessrelation to bring them to the desired form. First we re- place | v (cid:48) (cid:105)(cid:104) v (cid:48) | with ˆ − | c (cid:48) (cid:105)(cid:104) c (cid:48) | in (A7b) and | c (cid:48) (cid:105)(cid:104) c (cid:48) | withˆ − | v (cid:48) (cid:105)(cid:104) v (cid:48) | in (A7e). The two terms containing ˆ can bereduced to − r j Im [ (cid:104) r | c (cid:48) (cid:105)(cid:104) c (cid:48) | ˆ r l | c (cid:105)(cid:104) c | r (cid:105) ] = 0 , (A9)leaving − eE c E cv Im [ (cid:104) r | c (cid:48) (cid:105)(cid:104) c (cid:48) | ˆ r j | v (cid:105)(cid:104) v | ˆ r l | c (cid:105)(cid:104) c | r (cid:105) ] − eE v E cv Im [ (cid:104) r | v (cid:48) (cid:105)(cid:104) v (cid:48) | ˆ r j | c (cid:105)(cid:104) c | ˆ r l | v (cid:105)(cid:104) v | r (cid:105) ] . (A10)Next we combine the second term in (A7c) with the firstin (A7f) using ˆ H = E v (cid:48) | v (cid:48) (cid:105)(cid:104) v (cid:48) | + E c (cid:48) | c (cid:48) (cid:105)(cid:104) c (cid:48) | , eE cv Im (cid:104) (cid:104) r | ˆ H ˆ r j | v (cid:105)(cid:104) v | ˆ r l | c (cid:105)(cid:104) c | r (cid:105) + (cid:104) r | ˆ H ˆ r j | c (cid:105)(cid:104) c | ˆ r l | v (cid:105)(cid:104) v | r (cid:105) (cid:105) − eE cv Im [ E c (cid:48) (cid:104) r | c (cid:48) (cid:105)(cid:104) c (cid:48) | ˆ r j | v (cid:105)(cid:104) v | ˆ r l | c (cid:105)(cid:104) c | r (cid:105) + E v (cid:48) (cid:104) r | v (cid:48) (cid:105)(cid:104) v (cid:48) | ˆ r j | c (cid:105)(cid:104) c | ˆ r l | v (cid:105)(cid:104) v | r (cid:105) ] . (A11)Writing ˆ H ˆ r j as ˆ r j ˆ H − i (cid:126) ˆ v j and then canceling two terms according to Eq. (A9), the first line becomes − e (cid:126) E cv Re [ (cid:104) r | ˆ v j | v (cid:105)(cid:104) v | ˆ r l | c (cid:105)(cid:104) c | r (cid:105) + (cid:104) r | ˆ v j | c (cid:105)(cid:104) c | ˆ r l | v (cid:105)(cid:104) v | r (cid:105) ] . (A12)0Collecting terms in Eq. (A6) for the cross-gap local AHC we find, after some cancellations, σ AHcg ,i ( r ) = − e E cv (cid:15) ijl Re [ (cid:104) r | ˆ v j | v (cid:105)(cid:104) v | ˆ r l | c (cid:105)(cid:104) c | r (cid:105) + (cid:104) r | ˆ v j | c (cid:105)(cid:104) c | ˆ r l | v (cid:105)(cid:104) v | r (cid:105) ]+ πe h (cid:15) ijl Im [ (cid:104) r | v (cid:48) (cid:105)(cid:104) v (cid:48) | ˆ r j | c (cid:105)(cid:104) c | ˆ r l | v (cid:105)(cid:104) v | r (cid:105) − (cid:104) r | c (cid:48) (cid:105)(cid:104) c (cid:48) | ˆ r j | v (cid:105)(cid:104) v | ˆ r l | c (cid:105)(cid:104) c | r (cid:105) ] . (A13)Comparing the first line with Eq. (16) for ∂ E ˆ P and using projection operators in the second line, we arrive at σ AHcg ( r ) = e (cid:104) r | ˆ v × ∂ E ˆ P | r (cid:105) + πe h Im (cid:104) (cid:104) r | ˆ P ˆ r ˆ Q × ˆ Q ˆ r ˆ P | r (cid:105) − (cid:104) r | ˆ Q ˆ r ˆ P × ˆ P ˆ r ˆ Q | r (cid:105) (cid:105) . (A14)As noted earlier, the last two terms in this expression are equal to one another, resulting in Eq. (23). [1] D. Vanderbilt and R. D. King-Smith, “Electric polariza-tion as a bulk quantity and its relation to surface charge,”Phys. Rev. B , 4442 (1993).[2] R. D. King-Smith and David Vanderbilt, “Theory of po-larization of crystalline solids,” Phys. Rev. B , 1651(1993).[3] N. Marzari and D. Vanderbilt, “Maximally localized gen-eralized Wannier functions for composite energy bands,”Phys. Rev. B , 12847 (1997).[4] L. D. Landau and E. M. Lifshitz, Electrodynamics of Con-tinuous Media , 2nd ed. (Pergamon Press, Oxford, 1984).[5] M. Fiebig, “Revival of the magnetoelectric effect,” J.Phys. D: Appl. Phys. , R123 (2005).[6] X.-L. Qi, T. L. Hughes, and S.-C. Zhang, “Topologicalfield theory of time-reversal invariant insulators,” Phys.Rev. B , 195424 (2008).[7] A. M. Essin, J. E. Moore, and D. Vanderbilt, “Magne-toelectric Polarizability and Axion Electrodynamics inCrystalline Insulators,” Phys. Rev. Lett. , 146805(2009).[8] S. Coh, D. Vanderbilt, A. Malashevich, and I. Souza,“Chern-Simons orbital magnetoelectric coupling ingeneric insulators,” Phys. Rev. B , 085108 (2011).[9] M. Sitte, A. Rosch, E. Altman, and L. Fritz, “Topolog-ical Insulators in Magnetic Fields: Quantum Hall Effectand Edge Channels with a Nonquantized θ Term,” Phys.Rev. Lett. , 126807 (2012).[10] M. Taherinejad and D. Vanderbilt, “Adiabatic Pumpingof Chern-Simons Axion Coupling,” Phys. Rev. Lett. ,096401 (2015).[11] T. Olsen, M. Taherinejad, D. Vanderbilt, and I. Souza,“Surface theorem for the Chern-Simons axion coupling,” Phys. Rev. B , 075137 (2017).[12] F. Wilczek, “Two applications of axion electrodynamics,”Phys. Rev. Lett. , 1799 (1987).[13] A. Malashevich, I. Souza, S. Coh, and D. Vanderbilt,“Theory of orbital magnetoelectric response,” New J.Phys. , 053032 (2010).[14] A. M. Essin, A. M. Turner, J. E. Moore, and D. Van-derbilt, “Orbital magnetoelectric coupling in band insu-lators,” Phys. Rev. B , 205104 (2010).[15] D. Xiao, M.-C. Chang, and Q. Niu, “Berry phase ef-fects on electronic properties,” Rev. Mod. Phys. , 1959(2010).[16] A. Marrazzo and R. Resta, “Locality of the anomalousHall conductivity,” Phys. Rev. B , 121114 (2017).[17] R. Bianco and R. Resta, “Mapping topological order incoordinate space,” Phys. Rev. B , 241106 (2011).[18] M. Kohmoto, B. I. Halperin, and Y.-S. Wu, “QuantizedHall effect in 3D periodic systems,” Physica B , 30(1993).[19] E. I. Blount, “Formalisms of Band Theory,” Solid StatePhys. , 305 (1962).[20] T. B. Boykin, “Current density and continuity in dis-cretized models,” Eur. Phys. J. , 1077 (2010).[21] F. D. M. Haldane, “Model for a quantum Hall effect with-out Landau levels: Condensed-matter realization of the“parity anomaly”,” Phys. Rev. Lett. , 2015 (1988).[22] F. Schindler, A. M. Cook, M. G. Vergniory, Z. Wang,S. S. P. Parkin, B. A. Bernevig, and T. Neupert, “Higher-Order Topological Insulators,” Sci. Adv.4