Kinetic magnetoelectric effect in topological insulators

The kinetic magnetoelectric effect is an orbital analogue of the Edelstein effect and offers an additional degree of freedom to control magnetisation via the charge current. Here we theoretically propose a gigantic kinetic magnetoelectric effect in topological insulators and interpret the results in terms of topological surface currents. We construct a theory of the kinetic magnetoelectric effect for a surface Hamiltonian of a topological insulator, and show that it well describes the results by direct numerical calculation. This kinetic magnetoelectric effect depends on the details of the surface, meaning that it cannot be defined as a bulk quantity. We propose that Chern insulators and Z_2 topological insulators can be a platform with a large kinetic magnetoelectric effect, compared to metals by 5 - 8 orders of magnitude, because the current flows only along the surface. We demonstrate the presence of said effect in a topological insulator, identifying Cu$_2$ZnSnSe$_4$ as a potential candidate.

In the KME, the electric field induces the magnetization, which may look similar to the magnetoelectric (ME) effect [27][28][29][30][31][32][33][34][35] . Nonetheless, in the spin Edelstein effect and KME, metallic systems are considered, and nonequilibrium electron distribution by the electric field is a key to generate magnetization. On the other hand, the ME effect originates from the change of the electronic band structure by the electric field, while the electron distribution is assumed to stay in equilibrium. Thus the ME effect is mainly considered in insulators, but can also exist in metals 36,37 . In accordance with this differences in mechanisms between the KME and the ME effects, their symmetry requirements are different. In the ME effect, inversion and time-reversal symmetries should be broken. On the other hand, in the Edelstein effect, the inversion symmetry should be broken, but the time-reversal symmetry can either be preserved or borken. It is seen in chiral systems [1][2][3]10,11,13,14 , and in polar systems 15 .
Magnetoelectric tensors may require careful consideration of boundaries. While orbital magnetization is independent on the boundary 38-43 , the orbital magnetization when an electric field is applied may not have such properties. The general orbital magnetoelectric response [32][33][34][35] depends on the boundary 41 . Therefore, it is important to study the effect of the boundary of the response of orbital magnetization.
In this paper, we investigate KME in topological insulators such as three-dimensional Chern insulators and Z 2 topological insulators (Z 2 -TIs) [44][45][46] in which the currents are localized on the surface. First, we calculate the KME in three-dimensional topological insulators with chiral crystal structure. Second, we derive the KME based on the surface Hamiltonian, and we show that this effect depends on surface states. Finally, we propose candidate materials for this effect and estimate the values of the KME. By comparing the results with the results in a chiral semiconductor tellurium 5 , we show that in topological insulators the orbital magnetization as a response to the current is much larger than metals by many orders of magnitude. Results Formulation for KME. We consider a crystal in a shape of a cylinder along the z-axis, and calculate its orbital magnetization along the z-axis generated by the current along the z-axis. Let c be the lattice constant along the z-axis. We introduce the velocity operator v as v = − i [r, H], where r is the position operator and H is the Hamiltonian. In the limit of the system length along the z-axis to be infinity, the orbital magnetization at zero temperature is × − e 2 dxdyψ † n,kz (x, y)(r × v) z ψ n,kz (x, y), (1) where −e is the electron charge, ψ n,kz (x, y) and E n are the nth occupied eigenstates and energy eigenvalues of H at the Bloch wavenumber k z , respectively. f (E) is the distribution function at the energy E, S is the cross section of the crystal along the xy-plane and N is the number of occupied states. We note that the position a e d c b z (c) Brillouin zone of our model with high-symmetry points. (d, e) Energy bands for the Hamiltonian H with parameters (d) t x = t y = m = b x = b y , t 3 = 0.1t y and t 4 = 0.15t y and (e) t x = t y = m = b x = b y , t 3 = 0.001t y and t 4 = 0.15t y . The symbols such as Γ and X represent high-symmetry points in the wavevector space. (f) One-dimensional model with periodic boundary condition in z direction. It has a rectangle shape of a size L x × L y within the xy plane. The dots are the sites, and th blue squares denote the unit cells, as specified in (a). In order to see the boundary effect on KME, the outermost layers on the xz surface has no chiral hoppings and those on the yz surface has chiral hoppings. (g, h) KME calculated with parameters (g) t x = t y = m = b x = b y , t 3 = 0.1t y and t 4 = 0.15t y and (h) t x = b x = t y , t y = b y = m, t 4 = 0.15t y , and t 3 = 0.001t y . The blue, orange and green curves represent the data for L x = 10a, L x = 20a, and L x = 30a, respectively, while L y is fixed as L y = 30a. operator r is unbounded and problematic if a system is infinite in some directions. Nonetheless, in the present case, the system is infinite only along the z direction, while (r × v) z does not involve z, which means Eq. (1) is well defined. The unit of the orbital magnetization is [A/m] in the SI unit.
Then, within the Boltzman approximation 4 , the applied electric field E z changes f (E) from f 0 (E) into f (E) = f 0 (E) + eτ Ez ∂f 0 (E) ∂kz in a linear order in E z , where τ is the relaxation time assumed to be constant and f 0 (E) is the Fermi distribution function f 0 (E) = (e β(E−µ) + 1) −1 , β = 1/k B T, k B is the Boltzman constant and µ is the chemical potential. Then the orbital magnetization is generated as This is the KME. In Chern insulators, because the TRS is broken, an orbital magnetization is nonzero in equilibrium, and a current leads to a change in an orbital magnetization due to the KME. On the other hand, in Z 2 -TIs, an orbital magnetization vanishes in equilibirium due to TRS, and the KME leads to appearance of an orbital magnetization. We note that in addition to the off-equilibrium electron distribution discussed above, the electric field also modifies the electronic states. In Z 2 -TIs, where TRS is preserved, this modulation of the electronic states does not lead to the orbital magnetization in the linear order in E z , because the TRS is preserved within this order. On the other hand, in topological insulators without TRS, such as Chern insulators, this may also contributes to the orbital magnetization. This mechanism is similar to the conventional ME effect in insulators, and has been discussed also in metals 36,37 This calculation method is different from that in the previous study on metals 1,2,4 , where bulk contribution in a system infinite along x and y directions are calculated. Model calculation on a Chern insulator. As an example of a topological insulator, we consider an orbital magnetization in a Chern insulator with a chiral crystal structure. For this purpose, we introduce a threedimensional tight-binding model of a layered Chern insulator, as shown in Fig. 1a, connected via right-handed interlayer chiral hoppings (Fig. 1b). Each layer forms a square lattice within the xy-plane, with a lattice constant a, and they are stacked along the z-axis with a spacing c. as shown in detail in Methods. The Brillouin zone and the band structure is shown in Figs. 1c-e. We set the Fermi energy in the energy gap.
We calculate KME in a one-dimensional quadrangular prism with xz and yz surfaces shown in Fig. 1f (see Methods), with its results in Figs. 1g and 1h with the interlayer hopping t 3 = 0.1t y and t 3 = 0.001t y , respectively, for several values of the system size, L x and L y , representing the lengths of the crystal in the x and y directions. Thus, the KME is affected by boundaries and system size, and this size dependence remains even when the system size is much larger than the penetration depth of topological surface states. Therefore, KME cannot be defined as a bulk quantity. In contrast, the orbital magnetization in equilibrium is shown to be a bulk quantity in crystals [38][39][40] in 2005, which is nontrivial, because the operator for the orbital magnetization involves the position operator r. Later, we give an interpretation on this characteristic size dependence. Surface theory of KME for a slab. In topological insulators such as Chern insulators, only the topological surface states can carry a current. Here we calculate the KME using an effective Hamiltonian for the crystal surface. Thereby, we can capture natures of KME through this surface theory. We consider slab systems, with its surfaces on y = y ± (y + > y − ). The slab is sufficiently long along the x and z directions and we impose periodic boundary conditions in these directions. To induce the orbital magnetization M KME z,slab , we apply an electric field E z in the z direction. Due to the interlayer chiral hoppings, the surface current acquires a nonzero zcomponent (Fig. 2a).
) be the surfacestate dispersion on the y = y − surface as shown in Fig. 2b and Fig. 2c. For simplicity, we assume C 2z symmetry of the system. Then the surface state dispersion on the y = y + surface is given by Here, we assume that the surface states are sharply localized at y = y ± , namely, we ignore finite-size effects due to a finite penetration depth. Then, we rewrite equation (see Supplementary Note 1 for details). We note that the Fermi surface depends on the surface termination, and so does the KME. We also confirm the surface dependence from numerical calculations as shown in Figs. 2d-i. This formula applies to any topological insulators such as Z 2 -TIs 44-46 (see Supplementary Note 2). Surface theory of KME for a cylinder. From this slab calculation, we calculate the KME for a cylinder geometry. We consider a current along the z direction in a one-dimensional quadrangular prism with xz and yz surfaces (surfaces I-IV in Fig. 3a) through its surface Hamiltonian. Let L x and L y denote the system sizes along the x and y directions, respectively. Because the KME is sensitive to differences in crystal surfaces, as shown in slab systems, we consider the individual surfaces separately.
In particular, in Chern insulators we can calculate the energy eigenstates for the whole system from those for the surface Hamiltonians. For simplicity, we assume twofold rotation symmetry C 2z of the system, which relates between I and III, and between II and IV. Then, only the surface I and II are independent. We write down the eigenequations for these surfaces as where H I and H II are the surface Hamiltonians for the surfaces I and II, respectively, and ψ I kxkz = u I kxkz (k x , k z )e ikxx e ikzz and ψ II kykz = u II kykz (k y , k z )e ikyy e ikzz are Bloch eigenstates on the surfaces I and II, respectively.
We can determine these eigenstates from four conditions, equality of the energy eigenvalues, current conservation at the corner ( Fig. 3b) [47][48][49] , periodic boundary condition on the crystal surface and the normalization condition (see Supplementary Note 3).
Thus we obtain a formula for KME in a onedimensional prism of a three-dimensional Chern insulator When v x and v y are almost independent of k z , we approximate equation (6): 1t y , t 4 = 0.15t y and µ = 0. The red and blue curves represent Fermi surfaces on the y = y − and y = y + planes, respectively. (h, i) KME for the slab models I and II, shown as red solid lines and as blue broken lines, respectively, with parameters (h) t x = t y = m = b x = b y , t 3 = 0.1t y and t 4 = 0.15t y and (i) where M I,KME z,slab and M II,KME z,slab represent the KME for a slab (equation (3)) with the surface I and that with the surface II, respectively. Thus, the KME of the one-dimensional system can be well approximated by equation (7) expressed in terms of that for the slabs along xz and along yz planes.
In general topological insulators, we can also derive KME in terms of a simple picture of a combined circuit, consisting of four surfaces I-IV with anisotropic transport coefficients. We obtain where j circ is the circulating current density within the xy plane around the prism per unit length along the z-direction. σ I,II ij is the electric conductivity tensor for the surfaces I and II (see Supplementary Note 4). On the other hand, we can also show M I,KME z,slab = In Chern insulators, by using σ I xx ∝ v x and σ II yy ∝ v y , we arrive at equation (7). Thus, we can calculate the KME from the surface electrical conductivity from equation (8), which depends on the aspect ratio L x /L y .
We numerically comfirm that the results of direct calculation by equation (2) and those for surface calculation by equation (7) agree well (Fig. 4a-c). When the interlayer hopping is large (Fig. 4c), they slightly deviate from each other. This is because we cannot ignore the k z dependence of v x and v y and they are out of the scope of the approximate expression (7).
Finite-size effect. In our approximation theory, we assumed that the surface current is localized at the outermost sites and ignored a finite penetration depth. In fact, we can fit well the data with various system sizes One-dimensional prism of a Chern insulator. (a) Schematic figure of the one-dimensional prism of a Chern insulator, extended along z-axis, with its size L x × L y within the xy plane. We call the four side surfaces I, II, III and IV. (b) Schematic figure for current conservation at the corners around the one-dimensional prism in (a). Here, j I x and j II y are current densities in the circumferential direction of the prism on the surfaces I and II, respectively, and we impose them to be equal. v I x and v II y are corresponding velocities of electrons, and u I xz and u II yz are amplitudes of the wavefunctions in the respective surfaces.
with a trial fitting function which includes a finite-size effect in equation (7) (see Methods and Supplementary Note 5). From these results, the finite-size effect is of the order 1/L in the leading order, coming from the finite penetration depth. When the system size is much larger than the penetration depth, the result is well described by the surface theory as shown in Fig. 4d. Materials. Topological insulators without inversion symmetry can be a good platform for obtaining large KME, because the current flows on the surface. Therefore, the closed loop created by the current is macroscopic and it efficiently induces the orbital magnetization. In contrast, in the conventional KME in metals, a bulk current generates microscopic current loops in the bulk, which leads to a much smaller effect than a surface current in topological insulators. One can also regard this set of current loops as a macroscopic current loop along the surface, but the current in this case is of a microscopic amount, determined by the current per bulk unit cell. Therefore, the resulting effect in bulk metals is much smaller than the KME in topological insulators, where the current along the surface is of a macroscopic size. Moreover, the surface states of topological materials are robust against perturbations caused by impurities.
Under the non-inversion-symmetry constraint, we cannot diagnose Z 2 -TIs easily because the Z 2 topological invariant is expressed in terms of k-space integrals. Our idea here is to use S 4 symmetry to diagnose Z 2 -TIs, where we only need to calculate wavefunctions at four momenta according to the symmetry-based indicator theories [50][51][52] . After searching in the topological material database 53 , we notice that Cu 2 ZnSnSe 4 54 with 82 and CdGeAs 2 with 122 are two ideal candidates of Z 2 -TIs with a direct gap for obtaining a large KME (see Supplementary Note 6 for details). In the following, we will use Cu 2 ZnSnSe 4 , which a b d c

Surface Hamiltonian
Tight-binding model FIG. 4: Kinetic magnetoelectric effect (KME) in one-dimensional systems. (a-c) KME calculated from two different methods; one is a direct calculation by equation (2) (solid lines) and the other is by a combination of calculation results for surfaces along xz and yz planes based on equation (7) The results are shown for L x = 10a (blue), L x = 20a (orange) and L x = 30a (green), while L y is fixed to be 30a. (d) Dependence of the KME on the aspect ratio within the xy plane with parameters t x = b x = 1.1t y , t y = b y = m, t 4 = 0.15t y , t 3 = 0.001t y and µ = 0. Blue points represent the result of direct calculation by equation (2) for various system sizes. With parameter values t x = b x = 1.1t y , t y = b y = m, t 4 = 0.15t y , t 3 = 0.001t y and µ = 0, we calculate the KME with various system sizes from equation (2). The system sizes are (L x , L y ) ∈ {10a, 12a, . . . , 38a} × {10a, 12a, . . . , 38a}. Red points represent the result of fitting the numerical results of equation (2) with the fitting function in equation (13), and its limit for L x , L y → ∞ is shown as the dashed line. The solid line represents the results of the surface theory in equation (7).
only has S 4 symmetry as shown in Fig. 5a, as an example to show the magnitude of the KME with different surfaces and different surface terminations.
Since the magnetoelectric tensor for the space group can obtain an orbital magnetization M KME 1 by adding an external electric field E 1 , through the surface currents both on the [001] surface and on the [010] surface thanks to the nonzero α 11 (see Supplementary Note 7 for details). In our discussion of the KME effect, we set the current direction to be along z axis. Therefore we will set the 1-axis in the above magnetoelectric tensor to be the z axis in our theory.
Figures 5b and c are the Brillouin zone and the band structure of Cu 2 ZnSnSe 4 with a gap, through the first-principle calculations whose details are explained in Methods. On the [001] surface, terminations with Cu-Sn layer (surface A) and with Se layer (surface B) have different surface energies and Fermi surfaces, as shown in Fig. 5d-g, which contribute to a magnetoelectric susceptibility of α A 11 = −1.804 × 10 8 s −1 Ω −1 ·τ and α B 11 = −2.565 × 10 9 s −1 Ω −1 ·τ , respectively. On the A surface, there is a single surface Dirac cone at Γ point, forming an electron-like Fermi surface. On the B surface, the Dirac cone at Γ point forms an almost zero Fermi surface, but two surface Dirac cones at twoX momenta form two hole-like Femi surfaces. Because the Fermi surfaces on the B surface are much larger than those on the A surface, the magnetoelectric susceptibility on the B surface is one order of magnitude larger than that on the A surface. Similar calculations on the [010] surface are in the Supplementary Note 8, and the result is α C 11 = −2.324 × 10 9 s −1 Ω −1 ·τ for the surface C. Let us compare the results with metallic materials in the bulk. For simplicity, we focus on the cases with the electric field E and the resulting magnetization M KME along the z direction. This kinetic magnetoeletric response is expressed as M KME z = α zz E z , and the conductivity is j z = σ zz E z . Thus the magnetization in response to the current is M KME z = (α zz /σ zz )j z . In the relaxation time approximation, both α zz and σ zz are proportional to the relaxatoin time τ . For simplicity we consider the system to be a cube with its size L × L × L. In the bulk metallic systems, the current is carried by the bulk states, and α zz ∝ L 0 , σ zz ∝ L 0 . On the other hand, in topological systems, only the surface conducts the current, and the conductivity σ zz scales as σ zz ∝ L −1 . On the other hand, we have shown α zz ∝ L 0 , which means that α zz is an intensive quantity. Thus, in topological insulators, the scaling of the KME as a response to the electric field is represented by the response coefficient α zz ∝ L 0 . Meanwhile, the response coefficient α zz /σ zz to the current is proportional to L. It means that as a response to the current, topological materials will generate a large amount of orbital magnetization as compared to metals.
On the other hand, in topological insulators, the induced orbital magnetization as a response to the current becomes huge compared with metals. To show this, we consider a system with surfaces having anisotropic transport coefficients with sheet resistance σ xz, and σ zz, . Then we can define an angle θ by σ xz, /σ zz, = tan θ, where θ describes an angle between the electric field along the z direction and the surface current density j surf . For example, for the [001] surface of Cu 2 ZnSnSe 4 , we get σ A 21, 11, = 2.0 × 10 10 s −1 Ω −1 ·τ , which yield tan θ A = −0.49 and tan θ B = −0.26 by identifying x = 2 and z = 1. Then the total current along the z direction is 4Lj z while the circulating current is j circ = j surf x = j surf z tan θ. Thus the magnetization response M KME z to the current density j z (= 4Lj surf z /L 2 ) is M KME z /j = (j circ /j surf z )(L/4) = (L/4) tan θ, which is proportional to the system size. It is the response coefficient α/σ shown in Table I. Thus for the macroscopic system size, the response M KME z /j z is also of the macroscopic size such as milimeters, and it is many order of magnitude larger than that in tellurium, where M KME z /j z is evaluated to be M KME z /j z = 1.85 × 10 −9 m. This scaling to the system size shows a prominent difference in the KME in topological insulators from similar effects. In bulk metals studied in previous works, the KME is always independent of the system size L. In topological insulators, the current induces both the spin and the orbital magnetizations, and we propose that only the orbital magnetization shows a different scaling behavior. As a comparison, we calculate the spin magnetization induced by the electric field in Cu 2 ZnSnSe 4 , with calculation details and spin textures presented in Supplementary Note 9, and obtain α A,spin 11 = 2.798×10 −2 m·s −1 Ω −1 ·τ /L y , and α B,spin 11 = −2.575m · s −1 Ω −1 ·τ /L y . Thus it inversely scales with the system size along the y direction. Thus suppose the system size of 1mm, they are of the order 10 1 − 10 3 s −1 Ω −1 ·τ , and the orbital magnetization from the KME, having 10 8 − 10 9 s −1 Ω −1 ·τ , is larger than the spin counterpart by 5 -8 orders of magnitude. Thus, while the spin and orbital magnetization behave similarly in bulk metals, they are quite different in topological insulators, which is the main point of the present paper.

Conclusion
In summary, we propose KME in topological insulators with chiral structure. This KME is carried by surface current due to the asymmetric crystal structure of the surface. Therefore, the KME is sensitive to surface terminations, and it cannot be defined as a bulk quantity. We derive a formula for the KME as a surface quantity using the surface Hamiltonian, and show that it fits with numerical results.
In theoretical treatments, atomic orbitals can classify the orbital magnetization into intraatom and interatom contributions. Some atomic orbitals such as p x ± ip y have orbital angular momentum, which leads to cor- responding intraatomic orbital magnetization. On the other hand, the hopping between atoms lead, to the interatomic orbital magnetization. In tight-binding models with atomic orbitals, they are separately calculated. In some papers 10,11 , the intraatomic orbital magnetization is studied, while the interatomic one is studied in other papers 1,2,15 In real mateirals, these two contributions are not separable, and in the ab initio calculation 5 , their sum is calculated. In this paper, we found that in topological materials, the interatomic contribution is much larger due to the macroscopic current loop. We show that the response to the current in topological insulators is much larger than in metals. In Table I, we show scaling behaviors of the spin and intraatom/interatom orbital magnetizations for metals and TIs. The KME is shown as responses to an electric field, α, and to a current, α/σ. In particular, in the response to a current, α/σ, it scales as L 1 only in the interatom orbital magnetization, while other entities scales with L 0 . This shows a particular feature of the interatom orbital magnetization generated by a current proposed in the present paper.
In this paper, we put two Z 2 TIs as candidates for the topological KME. In addition to Chern insulators and Z 2 TIs, other classes of topological insulators such as various classes of topological crystalline insulators, and topolgical semimetals will also show the topological KME. In topological semimetals, topological surface states coexist with bulk metallic states, but the contribution from the former overwhelms that from the latter, and the topological KME is expected. Methods Details of the model Hamiltonian. We consider a Chern insulating system with a chiral crystal structure. The model is composed of infinite layers of the twodimensional Wilson-Dirac model 59,60 . The lattice sites are expressed by (i, j, l), with i, j, l being integers, specifying the x, y and z-coordinates. At each lattice site, we consider two orbitals 1 and 2. Let c i,j,l,σ denote the annihilation operator of electrons at the (i, j, l)-site with orbital σ (= 1, 2), and we write c i,j,l = (c i,j,l,1 , c i,j,l,2 ) T . t 4 ) , where H WD is an in-plane Wilson-Dirac Hamiltonian, and H interlayer is an interlayer Hamiltonian representing a structure similar to right handed solenoids. The in-plane Wilson-Dirac Hamiltonian is

The model Hamiltonian is
where H.c. stands for Hermitian conjugate of the preceding terms, † represents Hermitian conjugate, and m, t x , t y , b x and b y are real parameters. This Hamiltonian H WD can be rewritten in the momentum space as where k is the Bloch wavenumber. An isotropic version of the two-dimensional Wilson-Dirac model with b ≡ b x = b y and t x = t y exhibits the Chern insulating phase when 0 < m/b < 2 and 2 < m/b < 4 60 . Next we add interlayer hoppings, including a direct hopping t 4 along the z-axis and a chiral hopping t 3 , where t 3 and t 4 are real parameters. To describe the chiral hopping t 3 , the lattice sites in the square lattice in each layer into groups of four sites, (2i−1, 2j −1), (2i−1, 2j), (2i, 2j −1) and (2i, 2j) where i and j are integers, and we introduce chiral hoppings between the groups on the neighboring layers. Then the total Hamiltonian for this model on a tetragonal lattice is given by where These hoppings in H interlayer form structures similar to right-handed solenoids. When H WD in the Chern insulator phase, even if H WD is perturbed by H interlayer , the system remains in the Chern insulator with the Chern number within the xy plane equal to −1 as long as t 3 and t 4 are small. In the main text, we are interested in the KME due to the topological surface states in the topological Chern insulating phase, in which the Fermi energy is in the energy gap.
Fitting function for KME. By taking into account the finite-size effect in equation (7), we give a fitting fuction. The finite penetration depth of the surface states will lead to O(1/L) correction to the KME, and that around the corner will lead to O(1/L 2 ) correction 38 . Thus, the fitting function is where w i (i = 1, 2, . . . , 7) are real constants. I: Scaling behaviors of inter-and intraatom orbital and spin magnetizations versus the system size L, as a response to the electric field represented by α, and to the current α/σ in metals and in topological insulators (TIs). The scaling of the conductivity σ is also shown.

Supplementary Note 1. Surface theory in a slab of a 3D Chern insulator
In this section, for a slab of a 3D Chern insulator we explain how to rewrite Eq. (2) into Eq. (3) in the main text. Let us set the Chern number along the xy plane to be C xy = −1.
The energy of the chiral surface states on the two surfaces y = y ± based only on an intralayer Hamiltonian within the xy plane, representing the Chern insulator, is where k x is the wavevector along the x-direction and v x is the velocity along the x-direction.
It means that the surface current is along the x-direction. When we include the interlayer Hamiltonian having chiral hoppings, the surface current on the surface will acquire a nonzero z-component as shown in Fig. 2a. Therefore, the energies of the chiral surface states of this model on the surface y = y ± are Here, the velocity v z along the z-direction is added to represent chiral nature of the hopping along the surface. This simple form means that the velocity of the chiral surface states is a constant value (v x , v z ). Nonetheless, this formula needs to be modified by the following reason. The energy should be a periodic function of k z with a period of the reciprocal lattice vector G z = 2π c along the z-direction: E(k x , k z ) = E(k x , k z + 2π c ), but Eq. (S1) does not satisfy this requirement and we cannot use Eq. (S1). On the other hand, since E(k x , k z ) represents the topological chiral surface states, it is not periodic along the k x direction.
For example, on the y = y − surface, by increasing k x the surface states come out of the bulk valence band, and finally go into the bulk conduction band, as shown in Fig. 2b. Let E − = E − (k x , k z ) be the surface-state dispersion on the y = y − surface. For simplicity, we assume C 2z symmetry of the system. Then the surface state dispersion on the y = y + surface is given by Let ψ +(−) (k x , k z ) denote a wave function on the y = y + (y = y − ) surface of the crystal.
In this calculation, we assume that the surface states are sharply localized at y = y ± , which means that we ignore finite-size effects due to a finite penetration depth. From Eq. (3), we evaluate the KME as wherek = (k x , k z ) and v x,± (k) = 1 ∂E ± (k) ∂kx . Therefore, by noting that L y = y + − y − and v x,+ (k) = −v x,− (k), we obtain the KME as Henceforth we omit the subscript "−" in E − . In the limit of zero temperature T → 0, the k z derivative of the Fermi-Dirac function is expressed as where E l (k z ) and E u (k z ) are the lower and upper bounds of the surface state energies at k z .
By performing the integral over E, we obtain This is Eq. (3) in the main text for Chern insulators, with sgn ∂E(k x ,kz) ∂k x = +1 for the chiral surface states. We note that the right-hand side of this equation is not identically zero, because the derivative in the integrand is not a total derivative. By noting that we can rewrite Eq. (S9) as Likewise, the longitudinal conductivity is written as (S19) Supplementary Note 3. Surface Theory in a one-dimensional prism of 3D Chern insulator In this section, for a one-dimensional prism we explain how to rewrite Eq. (2) into Eq. (6).
We can determine the eigenstates on the surfaces I and II, u I kxkz and u II kxkz from four conditions. The first condition is that the energy eigenvalues are equal: E I kxkz = E II kykz . The second one is the current conservation at the corner 1-3 : where The third one is the periodic boundary condition on the crystal surface: where L x and L y are the lengths of the crystal in the x and y directions. This integer n can be regarded as a quantum number labeling the eigenstates, and we write their eigenvalues as E(n, k z )(= E I kxkz = E II kykz ). The last one is the normalization condition: From these conditions, u I kxkz and u II kykz are (S25) Then we obtain KME from Eq. (3) and Eq. (S25): When L x and L y are large, k x , k y and E(n, k z ) are regarded as continuous variables. By rewriting the summation over n to an integral over E and by using Eq. (S23), we obtain the KME as where we used From Eq. (S23), let k I0 x (k z , E) and k II0 y (k z , E) be the values of k x and k y for a given value of k z and E. Therefore they satisfy Therefore, we finally get: Therefore, we rewrite Eq. (S27) as where v I x = 1 ∂E I ∂kx and v II y = 1 ∂E II ∂ky and k I x (k z , E) and k II y (k z , E) are functions obtained from E = E I kxkz and E = E II kykz , respectively. This is the formula for KME in a one-dimensional prism of a three-dimensional Chern insulator.

Supplementary Note 4. KME and electric conductivity tensors
We can also derive KME in a Chern insulating system by using electric conductivity tensors on surfaces. On the surface I along the xz plane, the surface current density along xz plane, j x and j z in response to the electric field (E x , E z ), is represented as where σ I ij is the electric conductivity tensor. Since we assumed that the system is invariant under twofold rotation with respect to the z-axis, the surface III has the same tranport property with the surface I, but with an opposite normal vector. Namely, it has We note that the electric field along the z-axis, E z , is common among the four surfaces.
Similarly, on the surface II along the yz plane, the surface current density along the yz plane, j y and j z in response the electric field (E y , E z ), is written as Similarly, on the surface IV one has has We note that the electric field satisfies and that the current around the crystal is conserved: In this section, we introduce our fitting result in detail. We calculate the KME with various system sizes from Eq. (2). By using these results, we determine the parameters in the fitting function (see Methods). With the hopping parameters t x = b x = 1.1t y , t y = b y = m, t 4 = 0.15t y , t 3 = 0.001t y and µ = 0, we get the results fitting as where K ≡ e 2 τ E z t y / 2 , if we take the limit of L x → ∞ and L y → ∞ with the ratio L x /L y kept constant, we get: On the other hand, the result of Eq. (7) based on the surface theory is fitted as: Apart from the coefficients of L x in the numerators, Eqs. (S44) and (S45) agree well. Since the effect of L x in the numerators of Egs. (S44) and (S45) is much smaller than that of L y , it is difficult to numerically obtain the coefficient of L x in the numerator with good accuracy.
Thus, we conclude that these results agree well as shown in Fig. 4d.
As an example, we calculate α 33 in a one-dimensional quadrangular prism of Cu 2 ZnSnSe 4 from Eq. (8). From S 4 symmetry in Cu 2 ZnSnSe 4 , we get σ I 11 = σ II 22 and σ I 13 = −σ II 23 . Hence, when the quadrangular prism preserves the S 4 symmetry, i.e. when L x = L y , we get α 33 = 0 from Eq. (8), in accordance with Eq. (S46). On the other hand, when L 1 = L 2 , α 33 become nonzero as opposed to Eq. (S46), because the S 4 symmetry is broken by the shape of the sample.
Thus, in general systems, if a component of the kinetic magnetoelectric tensor α ij is nonzero from symmetry analysis, it means that the corresponding KME appears regardless of the crystal shape, reflecting the symmetry of the crystal. On the other hand, when α ij = 0, the KME does not occur when the crystal shape preserve the point-group symmetry of the crystal, but otherwise the KME may occur depending on the crystal shape. Compared with these values for the spin Edelstein effect, the orbital counterpart, i.e. the kinetic magnetoelectric effect, has α A 11, /σ A 11 = 0.25 and α B 11 /σ B 11, = 0.13, which are much larger by 6 orders of magnitude.