Half-quantized helical hinge currents in axion insulators

ABSTRACT Fractional quantization can emerge in noncorrelated systems due to the parity anomaly, while its condensed matter realization is a challenging problem. We propose that in axion insulators (AIs), parity anomaly manifests a unique fractional boundary excitation: the half-quantized helical hinge currents. These helical hinge currents microscopically originate from the lateral Goos-Hänchen (GH) shift of massless side-surface Dirac electrons that are totally reflected from the hinges. Meanwhile, due to the presence of the massive top and bottom surfaces of the AI, the helical current induced by the GH shift is half-quantized. The semiclassical wave packet analysis uncovers that the hinge current has a topological origin and its half quantization is robust to parameter variations. Lastly, we propose an experimentally feasible six-terminal device to identify the half-quantized hinge channels by measuring the nonreciprocal conductances. Our results advance the realization of the half-quantization and topological magnetoelectric responses in AIs.

To resolve the above puzzles, the microscopic mechanism of the boundary excitation in the AI should be clarified. So far, theoretical studies in the AI are mainly based on topological field theory, which cannot provide a definitive answer. Fortunately, the semiclassical analysis can give intuitive descriptions of microscopic processes in topological systems such as the adiabatic charge pumping, etc [34][35][36][37][38][39][40]. Specifically for the AI, electrons on the metallic side surfaces undergo multiple scattering between the insulating top and bottom surfaces. Therefore, it is highly desirable to illustrate the picture of these scattering processes based on semiclassical wave packet dynamics. This will further uncover the mechanism of boundary excitations in AIs and helps to provide a feasible proposal for the transport experiment to identify the AI phase.
In this work, we find that half-quantized helical hinge currents exist in the AI and manifest as its fingerprint. Based on the semiclassical wave packet dynamics, we establish the microscopic picture of these hinge currents. Similar to the Goos-Hänchen (GH) shift [41][42][43][44][45] of the light beam, the massless Dirac electrons on the side surfaces of the AI also undergo a lateral GH shift when reflected from the massive top or bottom surface (see Fig. 1). Moreover, due to the time-reversal (T ) symmetry-broken top and bottom surfaces, the GH shift favors a specific direction on the hinge. Consequently, helical GH shift currents are accumulated on the hinges of the AI when side-surface electrons bounce back and forth between the gapped top and bottom surfaces [see Fig. 1(a)], which is distinct from the chiral net currents on the edge of the Chern insulator (CI) [see Fig. 1(b)]. Interestingly, we find that the differential GH shift current δI GH is exactly half-quantized with respect to the differential Fermi energy δE F , i.e. δI GH = e/2hδE F . We demonstrate that the GH shift current originates from the non-vanishing Berry curvature during the scattering process, and its half-quantization is robust to the variation of parameters. Then, we numerically verify the half-quantized helical hinge currents through a 3D lattice model. Finally, we propose a six-terminal device to identify the half-quantized hinge currents in the AI by measuring the nonreciprocal conductances.
Model Hamiltonian and GH shift.-We model the AI by a 3D magnetic TI, of which the top and bottom surfaces are gapped oppositely while the side surfaces remain gapless. Such a model is in accordance with the experimental setups utilized to realize the AI state in magnetic TI heterostructure or antiferromagnetic TI MnBi 2 Te 4 [16][17][18]21]. According to the bulk-boundary correspondence theorem, the top/bottom surface and the side surface of the AI can be described by the 2D effective Hamiltonians: (1) Here, v F denotes the Fermi velocity, U is the gate potential, and m is the mass term induced by the magnetization. In the following discussions, we unfold the topside-bottom surfaces of the AI [see Fig. 1 We use the probability flux method [43,44] to calculate the GH shift on the hinge of the AI. As shown in Fig. 2(b), we place the massless side surface in x 0 and the massive top/bottom surface in x > 0. Firstly, we solve the scattering problem by matching ψ(r) = ψ in + rψ re (x 0) and ψ eva (r) = e −κx+ikyy ψ(0, 0) (x > 0) at x = 0, where ψ(r) is the superposition of the incident plane wave ψ in (r) = e ikxx+ikyy [e −i α 2 , e i α 2 ] T / √ 2 and reflected plane wave ψ re (r) = e −ikxx+ikyy [e −i π−α 2 , e i π−α 2 ] T / √ 2 in the massless region, and ψ eva (r) is the evanescent wave in the massive region. α = arctan ky kx denotes the incident angle [see Fig. 2 (a)] and r = (x, y). Here we only consider the total reflection process, which captures the physics within the magnetization gap. The reflection coefficient r = e iφr is given in [46]. Then, imagine that the incident and reflected waves have finite width as depicted in Fig. 2(b), the GH shift can be obtained under the probability flux conservation constraint. As labeled in Fig. 2(b), J int (J eva ) represents the flux carried by the interference wave (evanescent wave), J GH is the flux induced by GH shift, while J d is the flux through the cross-section colored blue with width d. The flux conservation condition imposes that J int + J eva = J GH + J d . Following the standard flux method analysis [43,44,46], we obtain In the upper panel of Fig. 2(c), we plot ∆ GH as a function of the incident angle α with fixed E, U , and m. For glancing incidence, i.e. α → ±π/2, ∆ GH diverges but with opposite sign for α = π/2 and α = −π/2. Such behavior indicates that at large incident angles, the GH shift has the same direction as the incident wave, thus does not show the chiral feature. However, ∆ GH peaks for vertical incidence, i.e. α → 0, showing clearly that the GH shift favors a specific direction and manifests the chiral feature. The chiral ∆ GH further implies that the net GH shift current accumulated on the hinge is chiral. We then decompose ∆ GH according to its contributions from the evanescent wave ∆ GH,eva and the interference wave ∆ GH,int [46]. As will be discussed later, ∆ GH,eva and ∆ GH,int introduce the components of the shift current with different decay laws. The lower panel of Fig. 2(c) shows that the chiral part of ∆ GH is mostly carried by ∆ GH,eva , which lies on the T -symmetry-broken top/bottom surface of the AI. In contrast, the non-chiral part, especially for glancing incidence α → ±π/2, is mostly carried by ∆ GH,int .
half-quantized GH shift current.-To give a microscopic picture of the chiral hinge current induced by the GH shift, we consider the electron in the scattering problem as a point particle, which bounces back and forth between double massive barriers as sketched in Fig. 2(a). Suppose the width between the opposite hinges is L x , thus the average time interval between two consecutive bounces is ∆τ = 2L x /v F cosα for −π/2 < α < π/2. When it bounces off the hinge, it undergoes a lateral GH shift ∆ GH . Such a lateral shift induces an anomalous velocity of electrons along the hinge as The total GH shift current induced by the v GH is obtained by counting the contributions of all filled electrons as where L y is the circumference of the side surface. K = (E + U )/ v F is the Fermi wave vector at energy E. Details for calculations can be found in [46]. By employing the stationary phase method [42,46], we find that the GH shift can be written as ∆ GH = −∂φ r /∂k y . Plugging this result into Eq. (4), we obtain the differential shift current δI GH = δE F [φ r (−π/2) − φ r (π/2)]e/2πh. In Fig. 2(d), we plot φ r as a function of α for different m. It is clear that φ r (−π/2) − φ r (π/2) = π and robust to the variation of m [46]. Therefore, δI GH is exactly half-quantized with respect to δE F as Moreover, δI GH can be decomposed according to its contributions from the evanescent wave and the interference wave, i.e. δI GH = δI GH,eva + δI GH,int . We prove in [46] that δI GH,int is power low (x − 1 2 ) decay from the hinge, while δI GH,eva decays exponentially from the hinge (e −x/λ ). Therefore, δI GH is different from the current carried by the topologically protected edge or hinge state, which decays exponentially on both sides of the hinge. The inset figure of Fig. 2(d) plots δI GH,eva and δI GH,int versus m. δI GH is mainly contributed from δI GH,eva when m is small. As m increases, the contribution from δI GH,int becomes dominant.
Topological origin of the GH shift current.-We provide a topological viewpoint of the half-quantized GH shift current in the frame of the adiabatic charge transport theory [34,36]. Here, we use Hamiltonian H(r) = v F (−iσ x ∂ x − iσ y ∂ y ) + m(x)σ z to describe the scattering process, where m(x) is a smooth function connecting the gapless and the gapped regions with m(x) → m for x 0 and m(x) → 0 for x 0. During the reflection process, the energy E and momentum k y are conserved. Therefore, the relation 2 v 2 F k 2 x + 2 v 2 F k 2 y +m 2 = E 2 holds. We take v F k x = −t as the virtual time and the scattering process is now described by the time-dependent 1D Hamiltonian H(k y , t) = −tσ x + v F k y σ y + m(t)σ z . The study of the GH shift is reduced to the study of an adiabatic charge transport problem in the y-direction. Define the Berry curvature Ω kyt = −2Im ∂u/∂k y |∂u/∂t , where |u(k y , t) represents the instantaneous eigenstate of H(k y , t). According to the adiabatic charge transport theory, v y (k y ) = ∂E(k y )/ ∂k y − Ω kyt is the velocity of the electron in the y-direction. The GH shift in the ydirection for a given k y is Combine it with Eq. (4), one obtains The integration of ∂E(k y )/ ∂k y vanishes because the band structure is symmetric with respect to k y . Note that Eq. (7) is nothing but the Berry phase along the boundary C of the integration manifold. As we have proved in [46], Γ(C) = −sgn(m)π, which is determined by the π Berry phase of massless Dirac electron and the sign of the massive barrier, thus is robust to local parameter changes.
Current density distribution on the hinges.-The halfquantized hinge current can be numerically visualized based on the 3D magnetic TI Hamiltonian. The model Hamiltonian is given by H = H 0 + H M , where H 0 = i=x,y,z Ak i τ x ⊗ σ i + (M 0 − Bk 2 )τ z ⊗ σ 0 is the nonmagnetic part and H M = M (r)τ 0 ⊗ σ z represents the magnetization term [47,48]. Here, σ i and τ i are the Pauli matrices acting on the spin and orbital spaces respectively. We use formula ImTr[ ∂H(k y ) ∂k y G r ky (E, r, r)]dk y (8) [31] to obtain the current density in the y-direction across the x-z plane with r = (x, z). Here, G r ky (E, r, r ) is the retarded Green's function for the momentum-sliced Hamiltonian H(k y ).
We first evaluate a semi-magnetic TI with a gapped top surface in y-z plane, which can be viewed as a "half" of an AI [see Fig. 3(a)] [27,28,49]. As shown in Fig. 3 Experimental characterization of the half-quantized helical hinge currents.-The half-quantized hinge current δI GH can be experimentally detected by a six-terminal device as illustrated in Fig. 4(a). Leads 1 and 3 (2 and 4) are contacted near the top (bottom) surface of the AI, while leads 5 and 6 are contacted to the ends of the sample. According to the Landauer-Büttiker formula, the transmission coefficient between terminals i and j is T ij = Tr[Γ i G r Γ j G a ] and the related differential conductance is G ij = e 2 /hT ij [50][51][52][53]. The measurement method of G ij is given in [46]. Γ i is the linewidth function of lead i and G r(a) is the retarded (advanced) Green's function of the AI sample.
To demonstrate the existence of the helical hinge channels, we calculate the nonreciprocal conductance G N ij = G ij −G ji between lead i and lead j. As shown in Fig. 4(c), in the AI state G N 13 = e 2 /2h and G N 24 = −e 2 /2h are halfquantized with opposite sign and G N 65 = 0. This implies that there exist two counterpropagating half-quantized hinge channels. Moreover, the spatial distribution of the local current J i→j (E) from site i to j [54] further reveals the helical signature of the transport hinge currents in the AI as shown in Fig. 4(b). Since the nonreciprocal conductance G N ij counts the asymmetric part of J i→j (E), we can deduce that the half-quantized conducting channel originates from the half-quantized chiral GH shift current δI GH on the hinge of the AI. The results are well consistent with those determined from GH shift (see Fig. 2) as well as the current density distribution (see Fig. 3). Therefore, these transport signals strongly confirm the reliability of the microscopic picture of half-quantized hinge currents proposed in Fig. 1. Since the half-quantized helical hinge currents serve as the fingerprint of the AI, our proposals also promote the experimental identification of the AIs.
Discussion.-The half-quantized helical hinge currents can be understood as a consequence of the quantized TME response in the AI [4]. As sketched in Fig. 4(d), a shift of the Fermi energy on the surfaces of the AI induces an interface electric field δE s = −∇δE F /e. According to the Maxwell equations modified by the axion term Lagrangian, δj = e 2 /2πh∇θ × δE s [4,9] with the spatially varying axion angle θ. The half-quantized hinge current δI = ±e/2hδE F is obtained by integrating δj over the corner.
Conclusion.-In conclusion, we found that the halfquantized helical hinge currents exist in AIs, and established their microscopic picture based on the GH shift currents. The half-quantization of the GH shift current has a topological origin and is robust to the variation of parameters. We numerically demonstrated that the halfquantized hinge channel is reflected by the half-quantized nonreciprocal conductances. Our studies deepen the perception of the boundary excitations of the AI and shed light on the detection of AIs through transport experiments. Sec7. Topological origin of the GH effect based on the adiabatic charge transport theory 7 Sec8. Derivation of the cross-section current density 9 Sec9. Derivation of the differential conductance and the local current density 10 Sec10. Chiral edge transport and its relation to half-quantized hinge channels in Chern insulators 13 Sec11. Experimental setups to measure the nonreciprocal conductances 14 Sec12. Model parameters in numerics 14

Sec1. INTRODUCTION FOR THE SUPPLEMENTARY MATERIALS
In this supplementary material, we give detailed explanations of the models and methods used in the main text, and give detailed derivations of the analytical results. In Sec2, we use the traditional stationary phase method to derive an analytical expression for the GH shift. In Sec3, we make necessary additions to the technical details of the probability flux method mentioned in the main text, and compare the evanescent wave and interference wave components of the GH shift as a function of E at different incidence angles. In Sec4, we derive the anomalous velocity induced by the GH shift and use a 2D lattice model to verify its modification to the band structure. In Sec5, we use the anomalous velocities induced by the GH shift to derive a half-quantized GH shift current in the AI. In Sec6, we discuss the power low decay of the interference wave part of the GH shift current, which serves as a unique feature of the half-quantized hinge current. In Sec7, we use the adiabatic charge transport theory to give a topological understanding of the halfquantized GH shift current. In Sec8, we derive the expression of the cross-section local current density. In Sec9, we use the non-equilibrium Green's function method to derive the multi-terminal differential conductance and the local current density. In Sec10, we additionally illustrate the relationship between chiral edge transport and half-quantized hinge channels in the CI, which is in stark contrast to the helical side surface transport in the AI phase. In Sec11, we propose a feasible experimental setup to elucidate the principle in measuring nonreciprocal conductances. Finally, in Sec12, we list the model parameters used in the numerical calculations.

Sec2. STATIONARY PHASE METHOD IN CALCULATING THE GH SHIFT
We use the stationary phase method to derive the expression of the GH shift ∆ GH [1][2][3]. We consider the scattering problem described by the 2D Dirac Hamiltonian used in the main text. The incident plane wave and the reflected plane wave are The scattering problem is solved by matching the boundary conditions of incident/reflected wave ψ(r) = ψ in + rψ re (x 0) and the penetrated evanescent wave ψ eva (r) = e −κx+ikyy ψ(0, 0) (x > 0), and obtain the reflection coefficient r = e iφr . We have and α = arctan ky kx . Therefore, φ r can be viewed as a function of k y or α for fixed E, U , and m.
The incident and reflected Gaussian wave packets are constructed as Note that the integration is only performed for k y because the eigen equation Hψ in/re g = Eψ in/re g constraints the number of free variables through k x = (E + U ) 2 − k 2 y and α = arctan ky kx are functions of k y . Expand α(k y ) and φ r (k y ) to the first order of (k y −k y ) aroundk y as α =ᾱ + (k y −k y )∂α/∂k y + O((k y −k y ) 2 ) and φ r = φ r + (k y −k y )∂φ r /∂k y + O((k y −k y ) 2 ). Then substitute them into Eq. (S3) and Eq. (S4) we have (accurate to the first order of k y ) for spin up (+) and spin down (−) components. From now on we replaceφ r ,k y , andᾱ with φ r , k y , and α. It is clear from Eq. (S5) and Eq. (S6) that the GH shift is the displacement of the wave packet center. GH shifts for spin up and spin down components are ∆ + GH = −∂φ r /∂k y − ∂α/∂k y and ∆ − GH = −∂φ r /∂k y + ∂α/∂k y . The total GH shift is the spin averaged wave packet center displacement ∆ GH = (∆ + GH + ∆ − GH )/2 = − ∂φr ∂ky . Substitute it into Eq. (S2) we have

Sec3. PROBABILITY FLUX METHOD IN CALCULATING THE GH SHIFT
The probability flux method is originally named "the energy flux method" in studying the GH shift of light beams based on the energy conservation condition [4,5]. Here, we treat the quantum wave system and replace the energy flux with the probability flux. Suppose the probability density of the incident/reflected beam is normalized to one, i. e. ψ † in (r)ψ in (r) = ψ † re (r)ψ re (r) = 1, then where ψ eva (r) = e −κx+ikyy ψ int (0, 0) for x > 0 and ψ int (r) = ψ in + rψ re for x 0. [See Fig. S1(a)]. The flux conservation condition implies that J int + J eva = J GH + J d . Following the standard flux method analysis, Therefore, the GH shift as a function of d can be obtained as The d dependence of ∆ GH (d) only comes from the last term in Eq . (S12), which oscillates with d and can be averaged out. Then we have It can be verified that Eq. (S7) and Eq. (S13) are equivalent. It is clear that the first term in Eq. (S13) comes from the evanescent wave and the second term comes from the interference wave. We plot ∆ GH,eva and ∆ GH,int as a function of E for different incident angles α in Fig. S1(b) and (c). We can see that for small incident angles (near vertical incidence), the total GH shift is mostly contributed from the evanescent wave while for large incident angles (near glancing incidence) the total GH shift is mostly contributed from the interference wave. Before deriving the half-quantized GH shift current in AIs and CIs, we investigate the anomalous velocity induced by the GH shift and its modification to the band structure. Consider a massless Dirac electron bounces back and forth between massive barriers with opposite mass as shown in Fig. S2(a). The average time interval between two consecutive bounces is ∆τ = Lx v F cosα . The anomalous velocity induced by the GH shift along the y direction is For electrons with incident angle α, the drift velocity (without the GH shift induced anomalous velocity) is v 0 = v F sinα. The total velocity vanishes at the band minimum v 0 + v GH = 0. We then have the condition for the band minimum The band minimum indicates the existence of the "8" shape trajectories as sketched in Fig. S2(a). We use the 2D tight binding Hamiltonian [6] to numerically investigate the band structure changes. We take y direction to be infinite and in the x direction the massless electrons with m = 0, µ = 0 sandwiched between two barriers with opposite m and µ = 0. The lattice Hamiltonian described by Eq. (S16) contains 4 valleys in total at (0, 0), (0, π), (π, 0), and (π, π). We investigate the two valleys at k y = 0 [(0, 0) and (0, π)]. Electrons from both valleys accumulate anomalous velocity when they bounce off the massive barrier, but with opposite direction as sketched in Fig. S2(a). From Fig. S2(b) we can see that band minimums for electrons from different valleys splits, which coincide with the analytical prediction based on Eq. (S15).

Sec5. DERIVATION OF THE HALF-QUANTIZED GH SHIFT CURRENT
We derive the half-quantized hinge current induced by the GH shift. Suppose the width between the barriers is L x , the average time between two consecutive bounces off one of the barriers is ∆τ = 2L x /v F |cosα| for −π < α π. (Here, we consider the time interval between two consecutive bounces of electrons off the same barrier, thus the distance traveled in the x direction is 2L x . Contributions from electrons with incident angle α and π − α are considered to be equivalent between successive bounces, i. e. ∆ GH (α) = ∆ GH (π − α), since the electron with incident angle α will be alternated to π − α after one reflection. Therefore, we take |cosα| to describe both cases.) The lateral GH shift ∆ GH induces an anomalous velocity of electrons near the barrier as v GH = ∆ GH /∆τ = ∆ GH v F |cosα|/2L x (−π < α π).
Some points should be noted for Eq. (S19) and Eq. (S20). In Eq. (S19), the integrand v F ∆ GH cosα is exactly the GH shift induced probability flux J GH in Eq. (S9). From the flux conservation perspective J GH = J eva + J int − J d , the GH shift current have both contributions from the evanescent wave and the interference wave, indicating that our derivations of the GH shift as well as the GH shift current are equivalent. In Eq. (S20), the integration over E is not necessarily performed from −∞ to E F , since we only derived the ∆ GH for the total reflection case within the gap. The only thing we are interested is the differential GH shift current with respect to E F , which contributes to the transport current. From Eq. (S7) ∆ GH = −∂φ r /∂k y and we have the differential form of I GH for E F lying in the gap of the top and bottom surfaces. From the plot in Fig. 2(d) in the main text, we immediately obtain the half-quantized GH shift current δI GH = e 2h δE F .

Sec6. POWER LAW DECAY OF THE INTERFERENCE WAVE PART OF THE GH SHIFT CURRENT
In the maintext and Sec3 we discussed the decomposition of the GH shift current according to the contributions from the evanescent wave part and the interference wave part. In this section, we emphasize that the interference wave induced GH shift current component decays from the edge following the power law, which is in stark contrast to the current carried by the topological edge or hinge state that decays exponentially. To see clearly, we now turn to the integrand in Eq. (S11). The y-component of the current as a function of x carried by the interference wave where k x = K cos α and K = (E F + U )/ v F is the Fermi wave vector. The total contribution of the current from all the α can be written as From Eq. (S2), if m → ∞, then φ r → 0 for most of the α. We omit the φ r for simplicity to study the asymptotic behavior of Eq. (S22). Then we have Here, J 0 (x) is the 0 th Bessel function. According to the asymptotic formula of the Bessel function one can see that when x → ∞ The result clearly shows that the interference wave part of the GH shift current maximizes at the boundary and decays to zero in a power law x − 1 2 when moving away from the boundary with the oscillation length as 1/2K. Note that the result in [7] is similar to ours, where they conclude that the edge current decays to zero in a power law x − 3 2 . The difference is that they counted the contributions from all states below the Fermi surface, while our result only focus on the states near the Fermi surface. Nevertheless, the results are consistent becasue the derivation of their result with respect to K d dK when x → ∞ is the same as ours.
Here we emphasize that the chiral current carried by the interference wave on the metallic side surface depends on the gapped, time-reversal symmetry breaking top/bottom surface. The power low decay of the current induced by the interference wave indicates that it cannot be generated by any kind of topologically protected edge or hinge state which decays exponentially. It also indicates that the half-quantized hinge current cannot exist by itself and should be combined with another one to form quantized side surface transport in the AI or the CI. These features are unique for the half-quantized hinge current.
Such a reflection process can also be understood as the 1D charge pumping problem in the y direction when we take vF kx = −t as the virtual time. The effective time-dependent Hamiltonian can be viewed as a Zeeman type H = B · σ with B = (−t, vF ky, m(t)) and the reflection process is simplified to an adiabatic spin procession for fixed ky and E. (c) Sketch of the charge pumping process in the y direction. The non-trivial Berry curvature Ω ky t induces a anomalous velocity v(ky) in the y direction. After the reflection, the total contribution from the adiabatic current gives rise to the GH shift.

Sec7. TOPOLOGICAL ORIGIN OF THE GH EFFECT BASED ON THE ADIABATIC CHARGE TRANSPORT THEORY
In this section, we give a semiclassical understanding of the topological origin of the GH effect based on the adiabatic charge transport theory [8][9][10][11]. Instead of Eq. (S1) that contains a sharp boundary between the massless and massive Dirac electron, we use Hamiltonian where m(x) is a smooth function connecting the gapless and the gapped regions with m(x) → m for x 0 and m(x) → 0 for x 0. Here we drop the U term in Eq. (S1) without affecting the conclusion. When the domain wall is large enough, the motion of the wave packet can be viewed semiclassically with specific momentum k and position r [see Fig. S4(a)]. The local Hamiltonian reads During the scattering process, the energy E and momentum in the y direction k y are unchanged. Therefore, the relation holds To map the scattering problem into a charge transfer problem, we take v F k x = −t as the virtual time. The local Hamiltonian in Eq. (S28) becomes time-dependent which describes the Zeeman coupling of a Pauli spinor to a time-dependent magnetic field B = (−t, v F k y , m(t)) [see Fig. S4(b)]. Then, the reflection process is reduced to an adiabatic spin procession. We denote the instantaneous eigen states of Eq. (S30) as |u ± (k y , t) , where ± denotes the spin up and spin down components of the spinor.
It is easy to solve the instantaneous eigenequation H(k y , t)|u ± (k y , t) = E ± |u ± (k y , t) and obtain |u ± (k y , t) = [E ± m, −t + i v F k y ] T / 2E(E ± m). Following the analysis in [11], apart from an unimportant overall phase factor up to the first order in the rate of the change in the Hamiltonian, the wave function is given by where n(n ) = ± represents the spin up or spin down components, and can further be viewed as the band index in a 1D (in the y direction) two-band model Eq. (S30). The average velocity for a given k y is found to the first order v n (k y )= ∂E n (k y )/ ∂k y − i n =n u n |∂H/∂k y |u n u n |∂u n /∂t Here, we used the relation u n |∂H/∂k y |u n = (E n − E n ) ∂u n /∂k y |u n and the identity n |u n u n | = 1. We only focus on the conduction band with n = + (the scattering process happens for electrons in the conduction band), thus from now on we omit the band index n. The GH shift in the y direction for a given k y is Fig. S4(c)]. From the maintext we show that δI GH /δE F = e h dky 2π ∆ GH . Then we have Since the band structure is symmetric with respect to k y , the integration of ∂E n (k y )/ ∂k y vanishes. Eq. (S36) is nothing but the Berry phase on the boundary C of the integration manifold. Intuitively, Γ(C) = π or Γ(C) = −π because on the boundary C we have T 2 + 2 v 2 F k 2 y = E 2 with m = 0. In such a case, C can be viewed as the Fermi surface of the massless Dirac cone, thus the Berry phase around C should be ±π. However, the sign of Γ(C) directly determines the direction of the chiral GH shift current according to Eq. (S36). To settle down this issue, the specific form of Ω kyt should be given.
Define k = t 2 + 2 v 2 F k 2 y and T = kcosθ, v F k y = ksinθ. The integral in Eq. (S36) can be performed in the polar coordinate system as where the Berry curvature Rewriting the spinor in the polar coordinate system as |u(k, θ) = [E + m, −ike −iθ ] T / 2E(E + m) and use the relation k 2 + m 2 = E 2 , we obtain the expression of Ω kθ after some derivations as Therefore, the Berry phase We conclude that the half-quantized chiral GH shift current δI GH /δE F = e 2h sgn(m) is protected by the π Berry phase of the massless Dirac electron while its direction is determined by the mass m of the massive barrier. Furthermore, it is quite easy to conclude that the result in Eq. (S40) is not affected by the potential U appeared in Eq. (S1), since it does not affect the π Berry phase and sgn(m). Nevertheless, the π Berry phase may also be influenced by the finite size gap or side surface random magnetization induced gap [12]. In these cases, the Berry phase is π(1 − δ E F ) where δ is the induced side surface gap. For very small δ, the Berry phase is approximately π and our analysis works well. Physically, the dependence of δI GH on m can be viewed as a consequence of the time-reversal symmetry breaking.
We make one more discussion on the relationship between the half-quantized current and the half-quantized charge pump as investigated in [11]. According to Eq. (2.6) in [11], the adiabatic charge pump is c n = −e T 0 dt BZ dq 2π Ω n qt , where T denotes the period of the cyclic pump and q denotes the momentum. Similarly, in our model the net charge pump during the reflection process is indicating that the net charge pump for the reflection on one of the massive barrier is exactly half-charge e 2 . When we take the other massive barrier to describe the side surface of the AI as depicted in Fig. S3(b) or the CI in Fig. S2(a), the two consecutive reflections on the two barriers make the charge pump process periodic. The total charge pumped along the y direction is 0 (1/2 − 1/2) in the AI and e (1/2 + 1/2) in the CI. The quantized e charge pump indicates the existence of the quantized chiral edge channel in CIs.

Sec8. DERIVATION OF THE CROSS-SECTION CURRENT DENSITY
The 3D TI Hamiltonian [13,14] used in calculating the cross-section current density and the conductances is where M (r) represents the spatially varying magnetization the couples to the TI through the Zeeman interaction.
Here, we take the system to be infinite in the y direction. The cross-section in the x − z plane is finite in x and y directions for the AI and the CI, and semi-infinite in the x direction for the semi-magnetic TI. Since the Hamiltonian is infinite in y direction for the three cases, k y is a good quantum number and the total Hamiltonian H can be decomposed into summations of the momentum-sliced Hamiltonians as where H(k y , r) is the momentum-sliced Hamiltonian with momentum k y where r = (x, z). We define the Green's function where |ψ n is the n th eigenstate of H(k y , r). Writing G r ky (E) in a real space form as G r ky (E, r, r ) = r|G r ky (E)|r . The velocity operator in the y direction for a given k y is v y (k y , r) = ∂H(k y , r)/ ∂k y . The local current density in the x − z plane can be expressed as where j y,ky (E, r) = − 1 π e ImTr[ ∂H(k y ) ∂k y G r ky (E, r, r)] is the local current density for a given k y . Eq. (S46) can be derived as follows j y,ky (E, r)= ev(k y , r)ρ(E, r) = n ev n (k y , r)δ(E − E ky,n ) (S47) In Eq. (S47), ρ(E, r) is the local density of states at energy E and position r, which can be expanded as ρ(E, r) = n δ(E − E ky,n ). In deriving Eq. (S49) from Eq. (S48), we used the relation 1 where P denotes the Cauchy principal value.

Sec9. DERIVATION OF THE DIFFERENTIAL CONDUCTANCE AND THE LOCAL CURRENT DENSITY
In this section, we derive the differential conductance and the local current density with the help of the nonequilibrium Green's function method [15][16][17]. As depicted in Fig. S5(c), we first consider the simplest case where the central region connects to an external lead. H C and H Lead are the Hamiltonians of the central region and the lead, H I is the coupling between the central region and the lead. We use a i to denote the annihilation operator in the central region at site i and c i to denote the annihilation operator in the lead at site i. The coupling Hamiltonian reads The particle leakage on the lead is Combine Eq. (S54) and Eq. (S55), we obtain the net current flowing into the lead as We denote the Green's function of the free lead (without coupling to the central region) as g, then use the Langreth theorem and Dyson equations we have Then we have We define self-energies and substitute it into Eq. (S59) as According to the fluctuation-dissipation theorem where A < = i(g r − g a ) is the spectral function and f Lead = 1/(e (E−µn)/k B T + 1) is the Fermi distribution function of the lead. We then have Define the linewidth function Γ Lead,i j = i(Σ r i j − Σ a i j ) and substitute it into Eq. (S61) Now we consider a multi-terminal system. We use Γ n and f n to denote the linewidth function and the Fermi distribution function of the n th lead. Define the total self-energy from all the leads According to the Keldysh formula we have The Dyson equations G r/a = g r/a + g r/a Σ r/a G r/a implies [g r/a ] −1 = [G r/a ] −1 + Σ r/a . Note that [g r ] −1 = [g a ] −1 = E − H Lead , therefore [G a ] −1 − [G r ] −1 = Σ r − Σ a = n −iΓ n . Then we have Combine Eq. (S64), Eq. (S66), and Eq. (S67), we obtain the current flowing into the m th terminal as At zero temperature, the Fermi distribution function becomes f n (E) = Θ(E F,n − E) = Θ(−eV n − E), where Θ(x) is the Heaviside step function and V n is the gate voltage applied to the terminal n. For small bias we take the approximation ∞ −∞ dE(f m − f n )(· · ·) ≈ (E F,m − E F,n )(· · ·) = e(V n − V m )(· · ·) and we have arrows in the upper panel of Fig. S5(d) and we take V 1 = V 2 to investigate the chiral or helical nature of the side surface current. In Fig. S5(e) we demonstrate that the nonreciprocal conductances on the top and bottom hinges of the CI have the same sign, contributing to totally quantized side surface transport. In Fig. S5(f) we illustrate that the quantized chiral conductance channel in the CI originates from the combination of the two half-quantized hinge channels (1/2+1/2).

Sec11. EXPERIMENTAL SETUPS TO MEASURE THE NONRECIPROCAL CONDUCTANCES
In this section, we illustrate the principles in measuring the nonreciprocal conductances in multi-terminal devices. We consider the nonreciprocal conductances between two external leads [such as lead 1 and lead 2 as shown in Fig. S5(g)]. To obtain G N 12 , the conductances G 12 and G 21 should be measured. G 12 is defined as G 12 = I 12 /V 1 , where V 1 is the applied gate voltage on lead 1 with all the other leads grounded as shown in Fig. S5(g), and I 12 is the current flowing into lead 2. Here, we only consider the differential conductances so that the I 12 is the differential current induced by the small bias V 1 . The measurement of G 21 = I 21 /V 2 is similar.
In experiments, fabricating the multi-terminal device in Fig. S5(g) (the bottom surface of the sample is grounded and all the other leads are connected near the top surface of the sample) may be easier than the six-terminal device shown in Fig. S5(f) [19]. Moreover, since the bottom surface of the sample is grounded in Fig. S5(f), the half-quantized hinge channel localized on the bottom hinge of the AI or CI is inactive, thus only the half-quantized hinge channel on the top hinge contributes to the nonreciprocal conductance. Therefore, the measurement result of G N 12 is sensitive to the quality of lead 1 and lead 2. To improve the experimental accuracy, we emphasize that the surface leads [lead 1 and lead 2 in Fig. S5(f)] should be thick enough (but do not touch the bottom surface) such that they couple to more conducting side surface channels. Besides, lead 1 and lead 2 should be close enough to ensure that nearly all hinge current can flow into the measuring lead.

Sec12. MODEL PARAMETERS IN NUMERICS
In Sec4, the parameters in the Hamiltonian Eq. (S16) are v F /a = 1. For the gapless region µ = 0.2 and m = 0, for the gapped region µ = 0 and m = 0.03. L x = 600.
In Sec10 and the main text, in calculating the cross-section local current density, we take the parameters of the Hamiltonian Eq. (S42) as A = 1, B = 0.6, M 0 = 1, E F = 0.4 and M = 0.6. In the semi-magnetic TI case, the magnetization term H M = 0.6 × τ 0 ⊗ σ z only couples to the top surface, while in the AI/CI case H M = ±0.6 × τ 0 ⊗ σ z couple to the top and bottom surfaces.
In calculating the local current distributions we take A = 1, B = 0.6, M 0 = 1, E F = 0.4 and M = 0.6. The system size (for both AI and CI) is 15 × 15 × 30. In plotting Fig. S5(d) and the Fig. 4(b) in the main text, the thickness is squeezed but does not affect the results in demonstrating the helical/chiral transport nature on the AI/CI side surface.