Magnetization energy current in the axial magnetic effect
MMagnetization energy current in the axial magnetic effect
Atsuo Shitade and Yasufumi Araki Institute for Molecular Science, Aichi 444-8585, Japan Advanced Science Research Center, Japan Atomic Energy Agency (JAEA), Tokai 319-1195, Japan (Dated: February 19, 2021)The axial magnetic effect (AME) is one of the anomalous transport phenomena in which theenergy current is induced by an axial magnetic field. Here we numerically study the AME for therelativistic Wilson fermion in the axial magnetic field and a twisted Dirac semimetal. The AMEcurrent density inside the bulk is nonzero, and particularly in the low energy regime for the formermodel, it is explained by the field-theoretical results without any fitting parameter. However, forboth models, the average AME current density vanishes owing to the surface contribution. The axialgauge field is regarded as the spatially modulated (effective) Zeeman field and induces the spatiallymodulated energy magnetization. The AME is attributed to the magnetization energy current andhence cannot be observed in transport experiments.
I. INTRODUCTION
The chiral anomaly has attracted much attention fromthe high energy and condensed matter communities.Chiral fermions have the chiral symmetry in the clas-sical action but not in the quantum-mechanical parti-tion function when parallel electromagnetic fields are ap-plied [1, 2]. Such systems are realized in quark-gluonplasmas in heavy-ion collision experiments [3–6] and ef-fectively in Dirac/Weyl semimetals [7].The chiral anomaly gives rise to various anomaloustransport phenomena. Among them is the chiral mag-netic effect (CME) in which the charge current is inducedby a magnetic field [8–12]. In the relativistic case, it isexpressed as j χ = q χµ χ B χ / π (cid:126) , (1a) j = j R + j L = q ( µ B + µ B ) / π (cid:126) , (1b) j = j R − j L = q ( µ B + µ B ) / π (cid:126) . (1c)Here, j R / L , µ R / L , and B R / L are the charge current den-sity, chemical potential, and magnetic field for right/left-handed fermions, respectively. q is the electric charge, µ ( µ ) = ( µ R ± µ L ) / B ( B ) = ( B R ± B L ) / µ cannot be realized [15, 16],as forbidden by the Bloch-Bohm theorem [17–19]. Byapplying parallel electric and magnetic fields, the chiralimbalance is generated away from equilibrium and resultsin the negative magnetoresistance via the CME [20, 21].In the condensed matter context, the phenomenon wasexperimentally observed in a noncentrosymmetric Weylsemimetal TaAs and its family [22–25].Another anomalous transport phenomenon is the chi-ral vortical effect (CVE) in which the charge current isinduced by the vorticity [10, 26–33]. In the relativisticcase, it is expressed as j χ = qχ ( µ χ + π T / ω / π (cid:126) , (2a) j = qµµ ω /π (cid:126) , (2b) j = q ( µ + µ + π T / ω / π (cid:126) . (2c)Here, T is the temperature, and ω is the vorticity. Sincethe rotating system of chiral fermions is in equilibrium,the CVE current is not expected to flow. However, theBloch-Bohm theorem cannot be applied to the CVE, be-cause the theorem is valid only in the thermodynamiclimit, while the system size is limited by the causal-ity [19, 27]. Recently, we found that the transport currentof the CVE vanishes regardless of the presence or absenceof µ [34]. In other words, the local charge current of theCVE is just the magnetization one that cannot be ob-served in transport experiments. We also demonstratedthat the anisotropic CVE can be observed in condensedmatter systems that belong to some chiral point groups.The axial magnetic effect (AME) is also an anomaloustransport phenomenon in which the energy current is in-duced by an axial magnetic field [35, 36]. In the relativis-tic case, the energy-momentum tensor is symmetric, andthe AME is reciprocal to the CVE in Eq. (2c). Hence, itis expressed as j ε = q ( µ + µ + π T / B / π (cid:126) . (3)The nonzero temperature part of Eq. (3), π T /
3, wasnumerically supported by the large-scale SU(2) latticegauge theory [36]. However, since the transport currentof the CVE vanishes [34], that of the AME should vanishas well. More generally, it is natural to expect that theenergy current should not flow in equilibrium, althoughthe Bloch-Bohm theorem for the energy current has notbeen proved yet.An axial gauge field can be effectively engineered inDirac/Weyl semimetals [13, 14, 37–39]. In magnetic Weylsemimetals, the Zeeman field from magnetic moments actas the axial gauge field, and a magnetic texture yields anaxial magnetic field [14, 37, 39]. Elastic strain can alsogenerate the axial gauge field in magnetic Weyl semimet-als [38] and nonmagnetic Dirac/Weyl semimetals [13, 14].Note that the axial gauge field in the latter case is spin de-pendent, because strain preserves the time-reversal sym- a r X i v : . [ c ond - m a t . m e s - h a ll ] F e b metry. Thus, magnetic or strained Dirac/Weyl semimet-als are ideal platforms to test the possibility for the aboveanomalous transport phenomena including the AME.In this paper, we numerically show that the AMEis canceled by the surface contribution. We considertwo different models; one is the relativistic Wilsonfermion [40] in an axial magnetic field, and the other is alattice model of a twisted Dirac semimetal Cd As [13].We calculate the charge and energy current densities im-posing the open boundary conditions. In the low energyregime for the former model, we find that the charge andenergy current densities inside the bulk are explained byEqs. (1b) and (3) without any fitting parameter. Forboth models, the current densities inside the bulk arenonzero, but the average ones vanish owing to the sur-face contributions. Such current distributions are ex-plained by the magnetization charge and energy currents,namely, the spatially modulated orbital and energy mag-netizations induced by the spatially modulated (effective)Zeeman field. In equilibrium, neither the CPME nor theAME can be observed in transport experiments. II. WILSON FERMION
The relativistic Wilson fermion is expressed as [40] H ( k ) = h ( k ) · ρ z σ + m ( k ) ρ x , (4)in which h x ( k ) = v sin k x a , h y ( k ) = v sin k y a , h z ( k ) = v sin k z a , and m ( k ) = m + r (3 − cos k x a − cos k y a − cos k z a )is the Wilson term. ρ and σ are the Pauli matrices for theparticle-hole and spin degrees of freedom, respectively,and form the Dirac representation of the Dirac matricesas α = ρ z σ and β = ρ x . Below, we choose v = 1 asthe energy unit and r = (2 / √ v . The lattice constantis a = 1.We introduce an axial gauge field in the Landau gauge,i.e., H ( x ) = − ( qva/ (cid:126) ) A y ( x ) σ y , (5)with A y ( x ) = B z x . Note that such a way of introducingthe axial gauge field is different from the conventionalway of introducing the vector gauge field using the Wil-son line. Since the Wilson term breaks the chiral symme-try, it is impossible to introduce the axial gauge field ina gauge-invariant way. We choose qB z a / (cid:126) = 1 × − .We diagonalize the total Hamiltonian [Eq. (4) +Eq. (5)] for the setup depicted in Fig. 1. The mass de-pends on x ; m = 0 in the massless region of N x − N m sites, while m = 8 r in the massive regions of N m sites.The axial gauge field is present in the massless and mas-sive regions continuously. We impose the open boundarycondition in the x axis and the periodic ones in the y and z axes. The number of the sites are N x = 300 and N m =50, leading to N x − N m = 200, and N y = N z = 192.Figure 2 shows the lowest positive eigenvalue for particlesand the highest negative one for antiparticles. In the low m =8 r m =8 rm =0 xA y =B z xN m N m N x -2 N m FIG. 1. Configuration of the mass m and axial gauge field A y . The massless region of N x − N m sites is sandwiched bythe shaded massive regions of N m sites. -2.5-2-1.5-1-0.5 0 0.5 1 1.5 2 2.5 -3 -2 -1 0 1 2 3 ε / v k y a or k z a k z a =0 k y a =0 FIG. 2. Lowest positive and highest negative eigenvalues ofthe total Hamiltonian [Eq. (4) + Eq. (5)] for the Wilsonfermion. Black solid and red dotted lines represent the eigen-values as functions of k y a ( k z a = 0) and k z a ( k y a = 0),respectively. energy regime, the eigenvalues are quadratic as functionsof k y a ( k z a = 0) and linear as functions of k z a ( k y a = 0),which is consistent with the chiral Landau levels in thepresence of the Wilson term.The charge and energy current densities are j z ( x ) = qa N y N z (cid:88) nk y k z v zn ( x, k y , k z ) f ( (cid:15) n ( k y , k z )) , (6a) j zε ( x ) = 1 a N y N z (cid:88) nk y k z v zn ( x, k y , k z ) (cid:15) n ( k y , k z ) × f ( (cid:15) n ( k y , k z )) , (6b)where v zn ( x, k y , k z ) = (cid:126) − (cid:104) u n ( k y , k z ) | ∂ k z H ( x, k y , k z ) | u n ( k y , k z ) (cid:105) .The distribution function is f ( (cid:15) ) = [ e ( (cid:15) − µ ) /T + 1] − for (cid:15) > f ( (cid:15) ) = − [ e ( | (cid:15) | + µ ) /T + 1] − for (cid:15) <
0. Theaverage current densities are obtained by j z ( ε )ave = 1 N x − N m (cid:88) x j z ( ε ) ( x ) . (7)In Fig. 3(a), the charge current density at the surfaceis opposite to that inside the bulk. In Figs. 3(b) and(c), the charge current density inside the bulk is com-pletely explained by the second term of Eq. (1b) withoutany fitting parameter, as far as the low energy regime( | µ/v | < . T /v < .
05) is concerned. However, theaverage charge current density vanishes for any chemicalpotential µ and temperature T [14]. The same is truefor the energy current density as shown in Fig. 4. Weconclude that Eqs. (1b) and (3) correctly describe theCPME and AME inside the bulk, respectively, but arecompletely canceled by the surface contributions.Let us comment on the previous numerical results onthe CPME in Ref. [14]. The authors considered a similarmodel of a Weyl semimetal using A y ( x ) = A y + B z x in Eq. (5). Hence, the spatial dependence of the chargecurrent density is different from ours, but the averageone vanishes as expected. They also showed the chargecurrent density at the center of the system but not relateit to Eq. (1b).We also comment on other previous numerical resultson the AME using the lattice gauge theory [35, 36].The authors used the overlap fermion that preservesthe modified chiral symmetry [41]. They obtained thenonzero temperature part of Eq. (3) only, because theopen boundary condition cannot be imposed in theirsetup. In our setup, the chiral symmetry is broken inthe high energy regime, which causes the deviation ofthe charge and energy current densities inside the bulkfrom the expected values. In the low energy regime, ourresults are consistent with the previous ones. III. LATTICE MODEL OF CD AS Cd As was predicted to be a Dirac semimetal by first-principles calculation [42] and soon later experimentallyconfirmed [43–46]. The effective model near the Γ pointis expressed as [13] H ( k ) = h ( k ) + h ( k ) · ρ σ z . (8)Here, h ( k ) = c + (2 c /c )(1 − cos k z c ) + (2 c /a )(2 − cos k x a − cos k y a ), h x ( k ) = ( v/a ) sin k x a , h y ( k ) =( v/a ) sin k y a , and h z ( k ) = m − (2 m /c )(1 − cos k z c ) − (2 m /a )(2 − cos k x a − cos k y a ). Two sets of the Paulimatrices, ρ and σ , represent the spin-orbit coupled states | / , ± / (cid:105) and | / , ± / (cid:105) . Since two sectors of σ z = ± σ z =+1 sector only and call the model 1 / As followingRef. [13]. If both sectors are considered, the charge andenergy currents below should read the spin current andspin-dependent energy current, respectively. The mate-rial parameters are c = − . c = 10 .
59 eV˚A , c = 11 . , m = 0 . m = 18 .
77 eV˚A , m = 13 . , and v = 0 .
889 eV˚A [47]. The latticeconstants are a = 12 .
67 ˚A and c = 25 .
48 ˚A.The Hamiltonian (8) is gapless at k = ∓ Q = [0 , , ∓ Q ]with Q satisfying m = (2 m /c )(1 − cos Qc ). Around these points, Eq. (8) is approximated as H ± ( δ k ) = ˜ c ∓ rv δk z + vδ k ⊥ · ρ ⊥ ± rvδk z ρ z , (9)where ˜ c = c + m c /m , v = vc /m , and r =(2 m /vc ) sin Qc . Using the above parameters, we ob-tain Qc = 0 . c = − . v = 0 .
502 eV˚A,and r = 1 . / As .Twisting 1 / As effectively generates a uniformaxial magnetic field [13]. Displacement is expressed as u ( x ) = θ ( z/L z ) x × e z , where θ is the twist angle, L z is the length in the z axis. Strain is then expressed as u xz ( x ⊥ ) = θy/ L z and u yz ( x ⊥ ) = − θx/ L z and coupledto electrons as H ( x ⊥ , k z ) = g ( v/a ) u ⊥ z ( x ⊥ ) · ρ ⊥ sin k z c = qv (cid:126) sin k z c sin Qc A ⊥ ( x ⊥ ) · ρ ⊥ . (10)Here, g is a dimensionless coupling constant, A ⊥ ( x ⊥ ) = B z [ − y, x ] /
2, and B z = − ( (cid:126) /qaL z ) gθ sin Qc . Around theWeyl points, the total Hamiltonian is approximated as H ± ( x ⊥ , δk z ) =˜ c ∓ rv δk z + v [ δ ˆ k ⊥ ∓ q A ⊥ ( x ⊥ ) / (cid:126) ] · ρ ⊥ ± rvδk z ρ z . (11)Thus, A and B effectively act as the axial gauge andmagnetic fields, respectively.We diagonalize the total Hamiltonian [Eq. (8) +Eq. (10)] imposing the open boundary conditions inthe x and y axes and the periodic one in the z axis.The number of the sites in the x , y , and z axes are N x = N y = 80 and N z = 192, respectively. We use gθ = 10, and the corresponding axial magnetic field is qB z a / (cid:126) = − . × − . Figure 5 shows the obtainedband structure. The color map represents the surfaceweight of the wave function | u n ( k z ) (cid:105) , ρ surf n ( k z ) = surf (cid:88) x ⊥ |(cid:104) x ⊥ | u n ( k z ) (cid:105)| . (12)Here, the summation is taken over the surface up to twosites. We find the bulk chiral Landau levels and the Fermiarc surface states.The charge and energy current densities are j z ( x ⊥ ) = qa cN z (cid:88) nk z v zn ( x ⊥ , k z ) f ( (cid:15) n ( k z )) , (13a) j zε ( x ⊥ ) = 1 a cN z (cid:88) nk z v zn ( x ⊥ , k z ) (cid:15) n ( k z ) f ( (cid:15) n ( k z )) , (13b)where v zn ( x ⊥ , k z ) = (cid:126) − (cid:104) u n ( k z ) | ∂ k z H ( x ⊥ , k z ) | u n ( k z ) (cid:105) ,and f ( (cid:15) ) = [ e ( (cid:15) − µ ) /T + 1] − is the Fermi distributionfunction. The average current densities are obtained by j z ( ε )ave = 1 N x N y (cid:88) x ⊥ j z ( ε ) ( x ⊥ ) . (14) -5x10 -5 -4x10 -5 -3x10 -5 -2x10 -5 -1x10 -5 -5 -150 -100 -50 0 50 100 150 (a) ¯ ha j z / q v x / a -5x10 -5 -4x10 -5 -3x10 -5 -2x10 -5 -1x10 -5 -5 -5 -5 -5 -5 -5 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 (a) (b) µ / v averagecenter -5x10 -7 -7 -6 -6 -6 -6 -6
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 (a) (b) (c) T / v averagecenter FIG. 3. (a) Spatial dependence of the charge current density for the Wilson fermion in the axial magnetic field at µ/v = T /v = 0 .
05 The shaded areas represent the massive regions. (b) The chemical potential dependence at
T /v = 0 .
05 and (c) thetemperature dependence at µ/v = 0 .
05 of the average charge current density (black cross) and that at the center (red square).The red dotted lines represent the second term of Eq. (1b) with qB z a / (cid:126) = 1 × − . -8x10 -6 -7x10 -6 -6x10 -6 -5x10 -6 -4x10 -6 -3x10 -6 -2x10 -6 -1x10 -6 -6 -150 -100 -50 0 50 100 150 (a) ¯ ha j ε z / v x / a -5x10 -6 -6 -5 -5 -5 -5 -5 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 (a) (b) µ / v averagecenter -2x10 -5 -1x10 -5 -5 -5 -5 -5 -5 -5 -5 -5 -5
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 (a) (b) (c) T / v averagecenter FIG. 4. (a) Spatial dependence of the energy current density for the Wilson fermion in the axial magnetic field at µ/v = T /v = 0 .
05 The shaded areas represent the massive regions. (b) The chemical potential dependence at
T /v = 0 .
05 and (c) thetemperature dependence at µ/v = 0 .
05 of the average energy current density (black cross) and that at the center (red square).The red dotted lines represent Eq. (3) with qB z a / (cid:126) = 1 × − . In Fig. 6(a), the charge current density at the surface isopposite to that inside the bulk. We show in Figs. 6(b)and (c) that the average charge current density vanishesfor any chemical potential µ and temperature T [13]. Thesame is true for the energy current density as shown inFig. 7. Both the CPME and AME are canceled by thesurface contributions and cannot be observed in trans-port experiments. IV. DISCUSSION AND SUMMARY
We attribute the current distributions obtained aboveto the magnetization currents. The axial “gauge” fieldin the above two models is a uniquely defined quantityand cannot be identified as a true one [16]. Rather, itis regarded as the (effective) Zeeman field as in Eqs. (5) and (10). In the presence of a spin-orbit coupling, theorbital magnetization is induced by the Zeeman field as M = χ os ( A /a ) [48]. We expect the energy magnetiza-tion M ε = χ es ( A /a ) as well. Now that the Zeeman fieldis spatially modulated, the orbital and energy magneti-zations are also spatially modulated and give rise to thecharge and energy current densities j ( ε ) = ∇ × M ( ε ) [49].However, such magnetization currents cannot be ob-served in transport experiments.To confirm the above scenario, we calculate the or-bital and energy magnetizations for twisted 1 / As as follows. Instead of the total Hamiltonian [Eq. (8) +Eq. (10)], we consider H ( k ) = H ( k ) + H ( x ⊥ , k z ) byfixing x ⊥ = x ⊥ . x ⊥ is just the parameter that controlsthe effective Zeeman field as A ⊥ /a = B z [ − y , x ] / a .We impose the periodic boundary conditions and calcu-late the orbital [50], heat, and energy magnetizations [51], -40-20 0 20 40 -3 -2 -1 0 1 2 3 ε [ m e V ] k z c ρ s u r f FIG. 5. Band structure of the total Hamiltonian [Eq. (8)+ Eq. (10)] for twisted 1 / As . The color map repre-sents the surface weight (12) up to two sites. The red dottedlines represent the band structure of the unperturbed Hamil-tonian (8) at k ⊥ = 0. M = q (cid:126) a cN x N y N z (cid:88) n k [ m n ( k ) f ( (cid:15) n ( k ))+ b n ( k ) f ( − ( (cid:15) n ( k ))] , (15a) M q = 1 (cid:126) a cN x N y N z (cid:88) n k { m n ( k )[ (cid:15) n ( k ) − µ ] f ( (cid:15) n ( k ))+ b n ( k ) f ( − ( (cid:15) n ( k )) } , (15b) M ε = M q + µ M /q. (15c)Here, the Berry curvature b n ( k ) and magnetic moment m n ( k ) are defined as b n ( k ) = i (cid:104) ∇ k u n ( k ) | × | ∇ k u n ( k ) (cid:105) , (16a) m n ( k ) =( − i/ (cid:104) ∇ k u n ( k ) |× [ (cid:15) n ( k ) − H ( k )] | ∇ k u n ( k ) (cid:105) , (16b)and f ( − ( (cid:15) ) and f ( − ( (cid:15) ) are f ( − ( (cid:15) ) = − (cid:90) ∞ (cid:15) dzf ( z ) , (17a) f ( − ( (cid:15) ) = − (cid:90) ∞ (cid:15) dz ( z − µ ) f ( z ) . (17b)The spin-orbit magnetic susceptibility χ os was alreadystudied in this way [52, 53]. In Fig. 8(a), the z compo-nent of the orbital magnetization is not so affected by theperturbative Zeeman field, because it mainly comes fromthe separation of Weyl points [54] determined by h z ( k )in Eq. (8). On the other hand, in Fig. 8(b), the x and y components are almost proportional to A ⊥ /a and showthe vortex structure. We also obtain the similar resultson the energy magnetization in Figs. 9(a, b). Here, wereinterpret x ⊥ as true coordinates and calculate the dis-cretized versions of j z ( ε ) = ∂ x M ( ε ) y − ∂ y M ( ε ) x . Such anadiabatic approximation is valid only inside the bulk. InFig. 8(c) and 9(c), we show the charge and energy cur-rent densities at the center. The magnetization currentsobtained here coincide with those obtained by imposingthe open boundary conditions in Sec. III. These resultsstrongly support our scenario of the magnetization cur-rents.To summarize, we have numerically investigated theCPME and AME using the relativistic Wilson fermionand a lattice model of a Dirac semimetal Cd As . Thecharge and energy current densities inside the bulk arecorrectly described by the previous results (the secondterm of Eq. (1b) and Eq. (3) in the relativistic case).However, the average current densities completely vanishowing to the surface contributions. The axial gauge fieldis regarded as the spatially modulated (effective) Zeemanfield and induces the spatially modulated orbital and en-ergy magnetizations. The CPME and AME currents arethe corresponding magnetization currents. Thus, it isimpossible to observe the CPME or AME in transportexperiments.What we called anomalous transport phenomena arenot transport phenomena in equilibrium. At the field-theoretical level, Eqs. (1b) is the covariant charge currentthat is not conserved. Conserved is the consistent chargecurrent obtained by adding the Bardeen-Zumino (BZ)polynomial [55]. As a consequence, the CME in the firstterm of Eq. (1b) does not occur [16]. On lattices, wealways consider the conserved charge current, and theCME does not occur [15]. The CVE current (2) is notcorrected by the BZ polynomial. However, it turned outto be the magnetization charge current [34]. The CPMEcurrent in the second term of Eq. (1b) and the AMEcurrent (3) are also the magnetization charge and energycurrents. These three currents do exist but cannot beobserved in transport experiments. ACKNOWLEDGMENTS
We thank M. Chernodub for asking a question thatmotivated this work when we submitted Ref. [34] and T.Hayata for valuable comments on our manuscript. Thiswork was supported by the Japan Society for the Pro-motion of Science KAKENHI (Grant No. JP18K13508).Y. A. is supported by the Leading Initiative for ExcellentYoung Researchers (LEADER). [1] S. L. Adler, Phys. Rev. , 2426 (1969). [2] J. S. Bell and R. Jackiw, Nuovo Cimento A , 47 (1969). (a) -40 -30 -20 -10 0 10 20 30 40 x / a -40-30-20-10 0 10 20 30 40 y / a -0.35-0.3-0.25-0.2-0.15-0.1-0.05 0 0.05 ¯ha j z / q [meV] -0.005 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045-20 -15 -10 -5 0 5 10 15 20 (a) (b) µ [meV] averagecenter -0.005 0 0.005 0.01 0.015 0.02 0.025 0 1 2 3 4 5 6 7 8 9 10 (a) (b) (c) T [meV] averagecenter FIG. 6. (a) Spatial dependence of the charge current density for twisted 1 / As at µ = 0 meV and T = 0 . T = 0 . µ = 0 meV of the average chargecurrent density (black cross) and that at the center (red square). (a) -40 -30 -20 -10 0 10 20 30 40 x / a -40-30-20-10 0 10 20 30 40 y / a -5 0 5 10 15 20 25 30 ¯ha j ε z [meV ] -1.8-1.6-1.4-1.2-1-0.8-0.6-0.4-0.2 0 0.2-20 -15 -10 -5 0 5 10 15 20 (a) (b) µ [meV] averagecenter -1.6-1.4-1.2-1-0.8-0.6-0.4-0.2 0 0.2 0 1 2 3 4 5 6 7 8 9 10 (a) (b) (c) T [meV] averagecenter FIG. 7. (a) Spatial dependence of the energy current density for twisted 1 / As at µ = 0 meV and T = 0 . T = 0 . µ = 0 meV of the average energycurrent density (black cross) and that at the center (red square). (a) -10 -5 0 5 10 x / a -10-5 0 5 10 y / a ¯haM z / q [meV] -10-5 0 5 10-10 -5 0 5 10 (a) (b) y / a x / a -1-0.8-0.6-0.4-0.2 0 0.2 0.4 0.6 0.8 1 ang l e / π ¯haM ⊥ / q [meV] (a) (b) (c) ¯ ha j z / q [ m e V ] T [meV] µ =-15 meV-10 meV5 meV0 meV5 meV10 meV15 meV FIG. 8. (a) z component and (b) the x and y components of the orbital magnetization for twisted 1 / As at µ = 0 meVand T = 0 . x ⊥ . (c) The temperature dependence of the charge current density at the center. Thesymbols correspond to the magnetization charge current calculated from the orbital magnetization, while the lines correspondto that obtained by imposing the open boundary conditions in Sec. III. (a) -10 -5 0 5 10 x / a -10-5 0 5 10 y / a -104.34-104.33-104.32-104.31-104.3-104.29-104.28 ¯haM ε z [meV ] -10-5 0 5 10-10 -5 0 5 10 (a) (b) y / a x / a -1-0.8-0.6-0.4-0.2 0 0.2 0.4 0.6 0.8 1 ang l e / π ¯haM ε⊥ [meV ] -1.6-1.4-1.2-1-0.8-0.6-0.4-0.2 0 0 1 2 3 4 5 6 7 8 9 10 (a) (b) (c) ¯ ha j ε z [ m e V ] T [meV] µ =-15 meV-10 meV5 meV0 meV5 meV10 meV15 meV FIG. 9. (a) z component and (b) the x and y components of the energy magnetization for twisted 1 / As at µ = 0 meVand T = 0 . x ⊥ . (c) The temperature dependence of the energy current density inside the bulk. Thesymbols correspond to the magnetization energy current calculated from the energy magnetization, while the lines correspondto that obtained by imposing the open boundary conditions in Sec. III.[3] I. Arsene et al. (BRAHMS Collaboration), Nucl. Phys.A , 1 (2005).[4] B. B. Back et al. (PHOBOS Collaboration), Nucl. Phys.A , 28 (2005).[5] J. Adams et al. (STAR Collaboration), Nucl. Phys. A , 102 (2005).[6] K. Adcox et al. (PHENIX Collaboration), Nucl. Phys. A , 184 (2005).[7] N. P. Armitage, E. J. Mele, and A. Vishwanath, Rev.Mod. Phys. , 015001 (2018).[8] A. Vilenkin, Phys. Rev. D , 3080 (1980).[9] D. Kharzeev, Phys. Lett. B , 260 (2006).[10] D. Kharzeev and A. Zhitnitsky, Nucl. Phys. A , 67(2007).[11] D. E. Kharzeev, L. D. McLerran, and H. J. Warringa,Nucl. Phys. A , 227 (2008).[12] K. Fukushima, D. E. Kharzeev, and H. J. Warringa,Phys. Rev. D , 074033 (2008).[13] D. I. Pikulin, A. Chen, and M. Franz, Phys. Rev. X ,041021 (2016).[14] A. G. Grushin, J. W. F. Venderbos, A. Vishwanath, andR. Ilan, Phys. Rev. X , 041046 (2016).[15] M. M. Vazifeh and M. Franz, Phys. Rev. Lett. ,027201 (2013).[16] K. Landsteiner, Acta Phys. Pol. B , 2617 (2016).[17] D. Bohm, Phys. Rev. , 502 (1949).[18] Y. Ohashi and T. Momoi, J. Phys. Soc. Jpn. , 3254(1996).[19] N. Yamamoto, Phys. Rev. D , 085011 (2015).[20] H. B. Nielsen and M. Ninomiya, Phys. Lett. B , 389(1983).[21] D. T. Son and B. Z. Spivak, Phys. Rev. B , 104412(2013).[22] X. Huang, L. Zhao, Y. Long, P. Wang, D. Chen, Z. Yang,H. Liang, M. Xue, H. Weng, Z. Fang, X. Dai, andG. Chen, Phys. Rev. X , 031023 (2015).[23] J. Du, H. Wang, Q. Chen, Q. Mao, R. Khan, B. Xu,Y. Zhou, Y. Zhang, J. Yang, B. Chen, C. Feng, andM. Fang, Sci. China Phys. Mech. Astron. , 657406(2016). [24] C.-L. Zhang, S.-Y. Xu, I. Belopolski, Z. Yuan, Z. Lin,B. Tong, G. Bian, N. Alidoust, C.-C. Lee, S.-M. Huang,T.-R. Chang, G. Chang, C.-H. Hsu, H.-T. Jeng, M. Ne-upane, D. S. Sanchez, H. Zheng, J. Wang, H. Lin,C. Zhang, H.-Z. Lu, S.-Q. Shen, T. Neupert, M. Z. Hasan,and S. Jia, Nat. Commun. , 10735 (2016).[25] Z. Wang, Y. Zheng, Z. Shen, Y. Lu, H. Fang, F. Sheng,Y. Zhou, X. Yang, Y. Li, C. Feng, and Z.-A. Xu, Phys.Rev. B , 121112(R) (2016).[26] A. Vilenkin, Phys. Rev. D , 1807 (1979).[27] A. Vilenkin, Phys. Rev. D , 2260 (1980).[28] J. Erdmenger, M. Haack, M. Kaminski, and A. Yarom,J. High Energy Phys. , 055.[29] N. Banerjee, J. Bhattacharya, S. Bhattacharyya,S. Dutta, R. Loganayagam, and P. Sur´owka, J. High En-ergy Phys. , 94.[30] D. T. Son and P. Sur´owka, Phys. Rev. Lett. , 191601(2009).[31] K. Landsteiner, E. Meg´ıas, and F. Pena-Benitez, Phys.Rev. Lett. , 021601 (2011).[32] K. Landsteiner, L. M. Eugenio Meg´ıas, and F. Pena-Benitez, J. High Energy Phys. , 121.[33] J.-Y. Chen, D. T. Son, M. A. Stephanov, H.-U. Yee, andY. Yin, Phys. Rev. Lett. , 182302 (2014).[34] A. Shitade, K. Mameda, and T. Hayata, Phys. Rev. B , 205201 (2020).[35] V. Braguta, M. N. Chernodub, K. Landsteiner, M. I.Polikarpov, and M. V. Ulybyshev, Phys. Rev. D ,071501(R) (2013).[36] P. V. Buividovich, J. Phys. Conf. Ser. , 012018(2015).[37] C.-X. Liu, P. Ye, and X.-L. Qi, Phys. Rev. B , 235306(2013); Phys. Rev. B , 119904(E) (2015).[38] A. Cortijo, Y. Ferreir´os, K. Landsteiner, and M. A. H.Vozmediano, Phys. Rev. Lett. , 177202 (2015).[39] Y. Araki, A. Yoshida, and K. Nomura, Phys. Rev. B ,115312 (2016).[40] K. G. Wilson, Phys. Rev. D , 2445 (1974).[41] H. Neuberger, Phys. Lett. B , 141 (1998).[42] Z. Wang, H. Weng, Q. Wu, X. Dai, and Z. Fang, Phys. Rev. B , 125427 (2013).[43] M. Neupane, S.-Y. Xu, R. Sankar, N. Alidoust, G. Bian,C. Liu, I. Belopolski, T.-R. Chang, H.-T. Jeng, H. Lin,A. Bansil, F. Chou, and M. Z. Hasan, Nat. Commun. ,3786 (2014).[44] Z. K. Liu, J. Jiang, B. Zhou, Z. J. Wang, Y. Zhang, H. M.Weng, D. Prabhakaran, S.-K. Mo, H. Peng, P. Dudin,T. Kim, M. Hoesch, Z. Fang, X. Dai, Z. X. Shen, D. L.Feng, Z. Hussain, and Y. L. Chen, Nat. Mater. , 677(2014).[45] S. Borisenko, Q. Gibson, D. Evtushinsky, V. Zabolotnyy,B. B¨uchner, and R. J. Cava, Phys. Rev. Lett. , 027603(2014).[46] S. Jeon, B. B. Zhou, A. Gyenis, B. E. Feldman, I. Kimchi,A. C. Potter, Q. D. Gibson, R. J. Cava, A. Vishwanath,and A. Yazdani, Nat. Mater. , 851 (2014). [47] J. Cano, B. Bradlyn, Z. Wang, M. Hirschberger, N. P.Ong, and B. A. Bernevig, Phys. Rev. B , 161306(R)(2017).[48] S. Murakami, Phys. Rev. Lett. , 236805 (2006).[49] K. Jensen, P. Kovtun, and A. Ritz, J. High Energy Phys. , 186.[50] J. Shi, G. Vignale, D. Xiao, and Q. Niu, Phys. Rev. Lett. , 197202 (2007).[51] A. Shitade, Prog. Theor. Exp. Phys. , 123I01 (2014).[52] R. Nakai and K. Nomura, Phys. Rev. B , 214434(2016).[53] Y. Ominato, S. Tatsumi, and K. Nomura, Phys. Rev. B , 085205 (2019).[54] M. Koshino and I. F. Hizbullah, Phys. Rev. B , 045201(2016); Phys. Rev. B , 209903(E) (2019).[55] W. A. Bardeen and B. Zumino, Nucl. Phys. B244