Magnetic susceptibility in three-dimensional nodal semimetals
aa r X i v : . [ c ond - m a t . m e s - h a ll ] O c t Magnetic susceptibility in three-dimensional nodal semimetals
Mikito Koshino and Intan Fatimah Hizbullah
Department of Physics, Tohoku University, Sendai 980-8578, Japan (Dated: October 9, 2015)We study the magnetic susceptibility in various three-dimensional gapless systems, including Diracand Weyl semimetals and a line-node semimetal. The susceptibility is decomposed into the orbitalterm, the spin term and also the spin-orbit cross term which is caused by the spin-orbit interaction.We show that the orbital susceptibility logarithmically diverges at the band touching energy in thepoint-node case, while it exhibits a stronger delta-function singularity in the line node case. Thespin-orbit cross term is shown to be paramagnetic in the electron side while diamagnetic in thehole side, in contrast with other two terms which are both even functions in Fermi energy. Thespin-orbit cross term in the nodal semimetal is found to be directly related to the chiral surfacecurrent induced by the topological surface modes.
PACS numbers: 75.20.-g, 73.20.At, 03.65.Vf
I. INTRODUCTION
Recently, a great deal of attention has been focusedon the three-dimensional (3D) materials having zero-gap electronic structure with nontrivial topology. TheDirac semimetal is a representative example of suchgapless materials, where four energy bands are touch-ing at an isolated point in the momentum space, un-der the protection of the time-reversal and the spatialinversion symmetry.
By breaking either symmetry,the four-fold degeneracy splits into a pair of doubly de-generate Weyl-nodes and then the system is called theWeyl semimetal.
The realizations of the Dirac andWeyl semimetal phases have been reported in the recentexperiments.
Theoretically, it is also possible to havethe band touching on a line in the momentum space, andsuch line-node band models and candidate systems wereproposed.
The response of the nodal semimetals tothe electromagnetic field attracts a great interest, andin particular, the exotic magnetoelectric effect relatedto its topological nature has been a topic of extensiveresearch.
In this paper, we study the magnetic susceptibility ofthe three-dimensional nodal semimetals and show thatit also exhibits various unusual properties. In two-dimenisional (2D) zero-gap linear band (e.g., graphene),it is known that the orbital susceptibility is expressed as adelta function of the Fermi energy, which diverges in thenegative (i.e., diamagnetic) direction at the band touch-ing point.
A similar calculation was also done for 3DDirac Hamiltonian, and there the susceptibility exhibits aweaker logarithmic singularity in the gapless limit. Herewe calculate the magnetic susceptibility in a wider vari-ety of three-dimensional gapless systems, including Diracand Weyl semimetals as well as the line-node semimetal.The calculation is not a simple extension of the previ-ous works in that we deal with the contribution from thespin degree of freedom. The 3D nodal band structureoften originates from a strong spin-orbit interaction, andin such a situation the magnetic susceptibility containsthe spin-orbit cross term in addition to the usual orbital susceptibility and the spin susceptibility. Here we con-sider a parameterized 4 × while it exhibits a stronger delta-functionsingularity in the line-node case. The spin susceptibility,which describes the Pauli and Van-Vleck paramagnetism,is found to be completely Fermi-energy independent nearthe degenerate points in the point-node case. The spin-orbit cross term is shown to be an anti-symmetric func-tion of the Fermi energy, which is paramagnetic in theelectron side and diamagnetic in the hole side. In thenodal semimetals, the spin-orbit cross term is found tobe directly related to the orbital magnetization inducedby the topological surface modes. Besides the continuummodel, we also consider the corresponding lattice modelto confirm the validity of the continuum approximation.The paper will be organized as follows. We introducethe continuum model and the lattice model for the point-node / line-node semimetals in Sec. II. The calculation ofthe magnetic susceptibility is presented in Sec. III wherethe point-node case and the line-node case are separatelyargued in Secs. IIIA and IIIB, respectively. A brief con-clusion is given in Sec. IV. II. MODEL HAMILTONIAN AND ENERGYSPECTRUMA. Point-node semimetal
We consider a 4 × H = vτ x ( σ · p ) + mτ z + bσ z = (cid:18) m + bσ z v σ · p v σ · p − m + bσ z (cid:19) , (1)where p = ( p x , p y , p z ) is the momentum vector, σ =( σ x , σ y , σ z ) is the Pauli matrices for the spin degree of FIG. 1. (Above) Energy spectra of ε sµ (0 , p y , p z ) for Dirac semimetal case ( m = b = 0), semiconducting case ( m = 1 , b = 0 . m = 0 . , b = 1) for the Hamiltonian Eq. (1). (Below) Corresponding Landau level spectra at themagnetic field ∆ B = 1. The thick red curve represents the zero-th Landau level, while the dashed curve non-existing zero-thlevel (see the text). freedom, ( τ x , τ y , τ z ) is the Pauli matrices for pseudospinscorresponding to, e.g., sublattices or atomic orbitals, and is the 2 × v , the mass parameter m , and the intrinsicZeeman field b which can arise in a magnetic materialbreaking time-reversal symmetry. By diagonalizing theHamiltonian, we obtain four energy bands, ε sµ ( p ) = s q m + b + v p + µ b p v p z + m , (2)where p = | p | , s = ± and µ = ± . The energy spectrum ε sµ (0 , p y , p z ) is plotted in Fig. 1 for (a) m = b = 0, (b) | m | > | b | and (c) | b | > | m | . The case of m = b = 0 cor-responds to the Dirac semimetal, where the spectrum iscomposed of a pair of degenerate linear bands touching at p = 0. The case | m | > | b | describes the semiconductingspectrum where the energy band is gapped in the range | E | < | m | − | b | . Finally, | b | > | m | represents the Weylsemimetal where the the middle two bands touch at apair of isolated point-nodes p = (0 , , ±√ b − m /v ).In the presence of the external magnetic field B in z direction, the Hamiltonian (1) becomes, H = vτ x ( σ · π ) + mτ z + (cid:18) b + 12 gµ B B (cid:19) σ z (3)where π = p + ( e/c ) A , A is vector potential to give B = ∇× A , g is g -factor and µ B = e ~ / (2 m c ) is the Bohrmagneton, and m is the bare electron mass. Here theexternal magnetic field B enters the Hamiltonian in twodifferent ways, one in the orbital part p +( e/c ) A , and theother in the spin Zeeman term gµ B Bσ z . In the originalDirac equation for a relativistic electron, the magneticfield appears only through A , and the spin Zeeman term gµ B B emerges out of it when being reduced to the low-energy quadratic Hamiltonian. In the present Dirac-likeequation in the solid state material, in contrast, we needto add gµ B B separately from A .The Landau levels can be found by using raising andlowering operators, π x + iπ y = ( √ ~ /l B ) a † and π x − iπ y = ( √ ~ /l B ) a , which operate on the usual Landau-level wave function φ n in such a way that aφ n = √ nφ n − and a † φ n = √ n + 1 φ n +1 . Here l B = p c ~ / ( eB ) is themagnetic length. For n ≥
1, the eigenfunction can bewritten as c φ n − |↑ , ↑i + c φ n |↓ , ↑i + c φ n − |↑ , ↓i + c φ n |↓ , ↓i , where | s, s ′ i represents | σ z = s i ⊗ | τ z = s ′ i . The Hamil-tonian matrix for the vector ( c , c , c , c ) then becomes H n ≥ = m + b vp z ∆ B √ n m − b ∆ B √ n − vp z vp z ∆ B √ n − m + b B √ n − vp z − m − b , (4)where ∆ B = √ ~ v/l B , (5)and b is actually b + gµ B B . This gives the four energylevels ε sµp z n = s q m + b + v p z + ∆ B n + µ b p v p z + m (6)where s = ± , µ = ± and n = 1 , , , · · · . The n = 0 is aspecial case where the eigenfunction is written as c φ n |↓ , ↑i + c φ n |↓ , ↓i . The Hamiltonian matrix for ( c , c ) be-comes H n =0 = (cid:18) m − b vp z vp z − m − b (cid:19) , (7)which gives the only two energy levels ε µp z = − b − µ p v p z + m , (8)with µ = ± and no s index.The Landau level spectrum is plotted against p z inFigs. 1 (d), (e) and (f) for m = b = 0, | m | > | b | and | b | > | m | , respectively. We note that the Landau lev-els of n ≥ s = ± branches. The n = 0 level ε µp z [Eq. (8); indicatedby the thick red curve] is not symmetric except for b = 0,because the inverted level − ε µp z does not actually exist(dashed curve). A set of the existing and non-existinglevels, ( ε µp z , − ε µp z ), can be regarded as the n = 0 sec-tor of Eq. (6). More precisely, Eq. (6) can be simplifiedto s | b + µ p v p z + m | at n = 0, and it is equal to eitherof ε µp z or − ε µp z . We define the index τ sµp z = ± ε sµp z = τ sµp z ε µp z . (9)Therefore, the complete set of the Landau level spec-trum including zero-th level is obtained by a single ex-pression of Eq. (6), where the index n of the sector( s, µ, p z ) runs from 0(1) when τ sµp z = +1( − FIG. 2. (Above) Energy spectrum of ε sµ (0 , p y , p z ) for theline-node semimetal described by Eq. (10) with b ′ = 1. (Be-low) Corresponding Landau level spectrum at the magneticfield ∆ B = 0 . B. Line-node semimetal
A line-node semimetal can be described by a similarHamiltonian but with a different constant matrix. Herewe consider a particular model described by , H = vτ x ( σ · p ) + b ′ τ z σ x = (cid:18) b ′ σ x v σ · p v σ · p − b ′ σ x (cid:19) , (10)The last term b ′ τ z σ x represents the spin Zeeman fieldparallel to x , of which direction is opposite between τ z =+ sector and τ z = − sector. The energy bands are givenby ε sµ ( p ) = s r v p x + h v q p y + p z + µb ′ i , (11)with indexes s = ± and µ = ± where the zero-energycontour becomes a circle given by p x = 0 and q p y + p z = b ′ /v . Figure 2(a) shows the band structure at p x = 0 FIG. 3. Energy spectra of ε sµ (0 , p y , p z ) of the lattice modelsfor (a) the Dirac semimetal [Eq. (15) with m = b = 0 and u = r = 1] and (b) the line-node semimetal [Eq. (17) with b ′ = 1 . u = r = 1] . plotted against ( p y , p z ), where the line node is indicatedby a solid circle. The spectrum is immediately gappedwhen p x goes away from 0.The Hamiltonian under the external magnetic field B in z axis becomes H = v ( σ · π ) τ x + b ′ τ z σ x + 12 gµ B Bσ z . (12)Similar to the point-node case, the Landau levels can beobtained by replacing π x ± iπ y with a † and a , respec-tively. The analytic expression is not available in thiscase and here we numerically obtain the Landau levelsby diagonalizing the Hamiltonian matrix with high in-dexes of n >
100 truncated. Figure 2(b) presents theobtained Landau levels plotted against p z , where we ne-glect the Zeeman term gµ B B for simplicity. We see thatthe zero-energy Landau level persists in the region of − b ′ < vp z < b ′ . It actually corresponds to the region in which the energy band Eq. (11) at fixed p z has pointnodes on p x p y plane, and each of the point nodes accom-modates the zero-th Landau level as in graphene.At p z = 0, in particular, the Hamiltonian Eq. (12) canbe transformed by a certain unitary transformation as H = vπ x σ x + ( vπ y − b ′ τ z ) σ y + 12 gµ B Bσ z , (13)which describes a pair of decoupled 2 × τ z = ± sectors. The Landau levels are obtained as ε sn = s q ∆ B n + b ( n = 1 , , · · · ) ,ε = − b, (14)where s = ± and b = gµ B B . C. Lattice model
In addition to the continuum model introduced in theprevious section, we also consider the lattice version ofthe point-node / line-node semimetals, to check the va-lidity of the continuum model in the susceptibility calcu-lation. For the point-node case, we introduce a Wilson-Dirac type cubic lattice model, H = uτ x ( σ x sin k x a + σ y sin k y a + σ z sin k z a )+ [ m + r (3 − cos k x a − cos k y a − cos k z a )] τ z + bσ z , (15)where vp i ( i = x, y, z ) of Eq. (1) is simply replaced by u sin k i a with the hopping energy u and the lattice spac-ing a . The extra term with r is introduced to gap out thepoint nodes other than the origin ( k x , k y , k z ) = (0 , , k i .Figure 3(a) shows the energy band ε sµ (0 , p y , p z ) of theDirac semimetal case given by Eq. (15) with m = b = 0and u = r = 1.In the real space, the Schr¨odinger equation becomes Eψ ( r ) = [( m + 3 r ) τ z + bσ z ] ψ ( r ) − r τ z X i = x,y,z [ ψ ( r + a i ) + ψ ( r − a i )] − iu τ x X i = x,y,z [ ψ ( r + a i ) − ψ ( r − a i )] , (16)where ψ ( r ) is the four-component wave function at thelattice point r , and a i is the unit lattice vector in i (= x, y, z ) direction. The external magnetic field can be in-cluded by appending the Peierls phase to the hoppingterms, and also adding the Zeeman energy gµ B B/ b .For the line-node case, the lattice model correspondingto Eq. (10) becomes H = uτ x ( σ x sin k x a + σ y sin k y a + σ z sin k z a )+ r (3 − cos k x a − cos k y a − cos k z a ) τ z + b ′ τ z σ x . (17)Energy spectrum ε sµ (0 , p y , p z ) with b ′ = 1 . u = r =1 is shown in Fig. 3(b). III. MAGNETIC SUSCEPTIBILITY
The magnetic susceptibility are obtained as in usualmanner by differentiating the thermodynamic potentialwith the external magnetic field B . Now B appears inthe orbital term p + ( e/c ) A and in the spin Zeeman term gµ B B in the Hamiltonian, and in the following, we for-mally treat these two magnetic fields for the orbital partand the spin part separately as B o and B s , respectively.The thermodynamic potential of the whole electronicsystem is written asΩ = − β X α ϕ ( ε α ) (18)where ϕ ( ε ) = ln[1 + exp − β ( ε − ζ )] , (19)where β = 1 / ( k B T ) with the termperature T , ε α is theenergy of the eigenstate α , and ζ is the chemical poten-tial. The magnetization is then written as a sum of thespin part and the orbital part, M = − ∂ Ω ∂B = − (cid:18) ∂ Ω ∂B s + ∂ Ω ∂B o (cid:19) ≡ M s + M o . (20)The magnetic susceptibility is then decomposed as χ = − ∂M∂B (cid:12)(cid:12)(cid:12)(cid:12) B =0 = − (cid:20) ∂ Ω ∂B + 2 ∂ Ω ∂B s ∂B o + ∂ Ω ∂B (cid:21) B s = B o =0 ≡ χ s + χ so + χ o . (21) A. Point-node semimetal
In the point-node case, the magnetic susceptibility canbe derived analytically using the Landau level spectrumEq. (6). The thermodynamic potential per unit volumeis explicitly written asΩ = − β πl B X sµp z ∞ X n = n ∗ ϕ [ ε sµp z ( x n )] , (22)where l B = p c ~ / ( eB o ) is the magnetic length and wedefined X sµp z = X s = ± X µ = ± Z ∞−∞ dp z π ~ , (23) ε sµp z ( x n ) = ε sµp z n = s q m + b + v p z + x n + µ b p v p z + m , (24) x n = n ∆ x, ∆ x = ∆ B = 2 ~ ev c B o , (25) and n ∗ = (cid:26) τ sµp z = +1 , τ sµp z = − . (26)The n ∗ specifies the beginning of the Landau level se-quence in the sector ( s, µ, p z ), and it is introduced toremove the non-existing zero-th level argued in Sec. II A.In Eq. (22), the summation over n can be transformedusing the Euler-Maclaurin type formula into a continuousintegral over variable x plus a series of the correctionterms. By using a notation F ( x ) = ϕ [ ε sµp z ( x )], we have ∞ X n = n ∗ F ( x n ) = 1∆ x "Z ∞ F ( x ) dx + τ sµp z F (0) ∆ x − F ′ (0) ∆ x
12 + O (∆ x ) , (27)where the second term in the bracket changes the signdepending on the starting index n ∗ .Noting that ∆ x ∝ B o , the thermodynamic potentialcan be expanded in powers of B o asΩ = Ω + λ B o + λ B + O ( B ) , (28)whereΩ = − β π ~ v X sµp z Z ∞ ϕ [ ε sµp z ( x )] dx, (29) λ = − β π ~ v ~ ev c X sµp z τ sµp z ϕ [ ε sµp z (0)] , (30) λ = − β π ~ v (cid:18) ~ ev c (cid:19) X sµp z ( − ∂∂x ϕ [ ε sµp z ( x )] (cid:12)(cid:12)(cid:12) x =0 . (31)Finally, the magnetization components are obtained as M s = − gµ B ∂ Ω ∂b , M o = − λ , (32)and the susceptibility components as χ s = − (cid:16) gµ B (cid:17) ∂ Ω ∂b , χ so = − gµ B ∂λ ∂b ,χ o = − λ , (33)where we used ∂/∂B s = ( gµ B / ∂/∂b .By using the explicit form of ε sµp z ( x ), we end up withthe formulas, χ s = 14 π ~ v (cid:16) gµ B (cid:17) X sµp z Z ∞ dx " − f ′ ( ε ) (cid:18) ∂ε∂b (cid:19) + f ( ε ) (cid:18) − ∂ ε∂b (cid:19) , (34) χ so = 14 π ~ v ~ ev c gµ B X sµp z f [ ε sµp z (0)] , (35) χ o = − π ~ v (cid:18) ~ ev c (cid:19) X sµp z f [ ε sµp z (0)] ε sµp z (0) , (36) (a) m = b = 0 (b) m = 1, b = 0.5 (c) m = 0.5, b = 1 Δχ s Δχ so m = 1, b = 0 ε F Δ χ s / χ s ∗ Δ χ s o / χ s o ∗ Δ χ o / χ o ∗ v p z Δχ s Δχ so Δχ o Δχ s Δχ so Δχ o ε F Δ χ s / χ s ∗ Δ χ s o / χ s o ∗ Δ χ o / χ o ∗ v p z ε F Δ χ s / χ s ∗ Δ χ s o / χ s o ∗ Δ χ o / χ o ∗ v p z Band structureBand structureBand structureBand structureBand structureBand structure Δχ o FIG. 4. Magnetic susceptibility in (a) the Dirac semimetal case ( m = b = 0) and (b) the semiconducting case ( m = 1 , b = 0 . m = 0 . , b = 1). In each column, the top panel is the 90 ◦ rotated band structure ε sµ (0 , , p z )and the second, third and fourth panels plot ∆ χ s , ∆ χ so and ∆ χ o against the Fermi energy. The dashed curves in (a) are for m = 1 , b = 0. where ε = ε sµp z ( x ), f ( ε ) = (cid:2) e β ( ε − ζ ) (cid:3) − is the Fermidistribution function, and we used the relation ϕ ′ ( ε ) = − βf ( ε ) and Eq. (9). The variable ζ in f ( ε ) is the chem-ical potential, and it will be denoted by ε F in the limitof T = 0.In the spin component χ s , the first term in the bracketis contributed from the states at the Fermi surface, andit corresponds to the conventional Pauli paramagnetism.The second term depends on all the occupied states belowthe Fermi energy, and it describes the Van-Vleck para-magnetism. The spin-orbit component χ so , Eq. (35), isformally proportional to the number of states below theFermi energy in the zero-th level ε sµp z (0), including bothexisting ( τ sµp z = +1) or non-existing ( −
1) branches. Forthe orbital part χ o , Eq. (36), the summation in the index s gives in the limit of T → X s = ± f [ ε sµp z (0)] ε sµp z (0) = θ ( | ε F | − | ε µp z | ) | ε µp z | , (37) where θ ( x ) is the step function returning 1(0) for x > < ε sµp z (0) = s | ε µp z | . Consid-ering the minus sign in front of Eq. (36), each sectorlabeled by ( µ, p z ) contributes to the diamagnetic suscep-tibility in the energy range −| ε µp z | < ε F < | ε µp z | , andits amplitude is inversely proportional to the width ofthe energy window. The delta function diamagnetism ingraphene is obtained in the limit of ε µp z → Now let us apply the above susceptibility formula toseveral different situations with specific ( m, b )’s. Whencalculating χ , we introduce a momentum cut-off in theintegral as | p z | < ε c /v and x < ε c with a sufficiently largeenergy ε c , to avoid the divergence of the integral and alsoto simulate the finite-sized Brillouin zone in real materi-als. As a result, the calculated susceptibility is shown toinclude a constant term (i.e., independent of the Fermienergy) which explicitly depends on the cut-off ε c . Inthe real material, however, such an offset term should bedetermined by the whole band structure beyond the de- FIG. 5. (a) Surface state bands in the Weyl semimetalwith m = 0 and a finite b . L and R indicate the left andright surface states, respectively. (b) Corresponding surfacecurrent on xy plane. scription of the present linear Hamiltonian, and it is notdeterminate in the continuum scheme. In the following,therefore, we argue about only the relative susceptibilityas a function of the Fermi energy, neglecting the offsetterm. The problem of the offset term will be addressedlater in this section using the lattice model.We first consider the Dirac semimetal case ( m = b =0). The susceptibility at T = 0 is explicitly calculated asfunctions of the Fermi energy ε F as∆ χ s = 0 , ∆ χ so = e vπ ~ c gµ B c ~ ev ε F , (38)∆ χ o = − e v π ~ c ln ε c | ε F | , where ∆ χ s etc. represents the relative susceptibility withthe constant term appropriately chosen. The three sus-ceptibility components are plotted in separate panels inFig. 4(a), together with the 90 ◦ -rotated band structurein the top panel. The orbital diamagnetism ∆ χ o loga-rithmically diverges in the diamagnetic direction at theband touching point. The spin term ∆ χ s becomes com-pletely Fermi-energy independent, where the energy de-pendences of Pauli term and Van-Vleck term completely cancel out. The spin-orbital term ∆ χ so is a simple linearfunction monotonically rising in increasing ε F .The vanishing ∆ χ s can be clearly understood by con-sidering the process to increase the spin Zeeman term b . The Hamiltonian Eq. (1) with m = 0 is decoupled to τ x = ± sectors as H ± = ± v [ σ x p x + σ y p y + σ z ( p z ∓ b/v )] , (39)which describes the Weyl semimetal composed of inde-pendent Weyl cones centered at p = (0 , , ± b/v ). Thespin field b just horizontally separates the Dirac cones inthe wave-space, so that the spin density never changes inthis process.On the other-hand, the finite ∆ χ so indicates that thefinite orbital magnetization M o is induced by increasing b , and this is explicitly calculated as M o = 12 ∆ χ so bgµ B / e π ~ vc ε F b. (40)The physical origin of the induced M o becomes clear byintroducing the surface to the system. When we con-sider the slab geometry which is finite in x -directionwhile infinite in y and z , we have a pair of surface local-ized bands which connect the two Weyl cones as shownin Fig. 5(a). The surface states then give the counterflows parallel to y directions on the opposite surfaces asillustrated in Fig. 5(b), and they contribute to a magne-tization in z direction. The velocity of the surface statealong y is equal to v . The area density of the surfaceelectrons existing between the zero energy to the Fermienergy is n surf = ( ε F /v )(2 b/v ) / (2 π ~ ) by considering thecorresponding momentum space. The total surface cur-rent density then becomes j = evn surf = eε F b/ (2 π ~ v ).The magnetization induced by the surface current, j/c ,exactly coincide with Eq. (40).The susceptibility calculation can be extended to thegeneral cases with a finite mass m . At b = 0, in paticular,we have∆ χ s = e vπ ~ c (cid:16) gµ B c ~ ev (cid:17) × m θ ( | ε F | − m ) ln p ε F − m + | ε F | m , (41)∆ χ so = e vπ ~ c gµ B c ~ ev × sgn( ε F ) θ ( | ε F | − m ) q ε F − m , (42)∆ χ o = − e v π ~ c × ln 2 ε c m | ε F | < m, ln 2 ε c p ε F − m + | ε F | | ε F | > m, (43)which are plotted as dashed curves in Fig. 4(a) ( m =1 , b = 0). In ∆ χ s , the cancelation of Pauli term andVan-Vleck term becomes incomplete, leaving a positivecomponent inside the band. We also see the log peak of∆ χ o is truncated, and the increase of ∆ χ so is interruptedby the energy gap. Now the system is in the trivial phasewithout the surface states so ∆ χ so is contributed by thebulk states.The susceptibility formula for the general ( m, b ) is com-plicated and the full expression is presented in AppendixA. Here we plot some representative cases in Fig. 4(b)and (c). When we introduce b smaller than m [Fig. 4(b)( m = 1 , b = 0 . m = 1 , b = 0 in Fig. 4(a). When b is larger than m , the system enters the Weyl semimetalphase [Fig. 4(c) ( m = 0 . , b = 1)], and we observe asimilar behavior to the Dirac semimetal case in the low-energy region reflecting the emergence of the point nodes;i.e., the vanishing ∆ χ s , the linear increase of ∆ χ so andlog divergence of ∆ χ o . Outside the linear band region,∆ χ s becomes negative (i.e., relatively diamagnetic com-pared to the zero energy) near the second band edges ε F = ± ( m + b ), and it increases again outside. We noticethat ∆ χ so is always an odd function of ε F while othertwo components, χ s and χ o , are even functions. This isclosely related to the fact that the electron-hole symme-try of the energy spectrum is broken solely by the zero-thLandau level when B o and B s coexist.In Fig. 4, the three susceptibility components are plot-ted in different characteristic scales, χ ∗ s = e vπ ~ c η , χ ∗ so = e vπ ~ c η, χ ∗ o = e vπ ~ c , (44)where η = gµ B c ~ ev ε , (45)and ε is the energy unit to measure m , b and ε F . Therelative magnitudes of three components are determinedby the dimensionless factor η , which is the ratio of thespin magnetic moment gµ B / µ = ~ ev / ( ε c ) associated withthe energy scale ε . In the typical velocity operator v ofthe order of 10 m/s, for instance, η becomes the orderof ε / (100meV). When the typical energy scales for theband structure such as ε F and m are much smaller than100 meV, χ o is dominant, χ so is in the middle range, and χ s is the smallest. In this case, the system should bediamagnetic in total.Finally, we calculate the susceptibility of the corre-sponding lattice model Eq. (15) to check the validityof the continuum approximation. Here the computa-tions are done numerically; we first obtain the eigenener-gies of the Hamiltonian at the several discrete points of( B orb , B spin ) which are small enough, and compute thethermodynamic potential Ω at each point. Then we takethe derivative of Ω using the differential approximation,and obtain the susceptibility using Eq. (21). In the cal-culation, we assume a finite temperature ( k B T = 0 .
1) tosmear out the discrete level structure at finite B orb .The results for the Dirac semimetal [ m = b = 0 and u = r = 1] are shown as solid curves in Fig. 6, where (a) χ s (b) χ so ε F Point-node (c) χ o Line-node -0.5 χ s / χ s ∗ χ s o / χ s o ∗ χ o / χ o ∗ FIG. 6. Susceptibility (a) χ s , (b) χ so and (c) χ o calculatedfor (a) the Dirac semimetal [Eq. (15) with m = b = 0 and u = r = 1] and (b) the line-node semimetal [Eq. (17) with b ′ = 1 . u = r = 1] . (a), (b) and (c) plot χ s , χ so and χ o , respectively. Theyare the absolute χ ’s and not the relative ∆ χ ’s. The ver-tical axis is scaled by Eq. (44) with v replaced with thecorresponding quantity, ua/ ~ . In the low-energy region | ε F | <
1, we actually see that the energy dependence ofeach χ well resembles the continuum counterpart in Fig.4(a), including the flat dependence of χ s , the linear be-havior χ so and the negative log peak of χ o . We see that χ s and χ o are even function, and both have the para-magnetic offsets compared to the continuum model. Incontrast, χ so is an odd function and it vanishes at ε F = 0.The results clearly show that the continuum calculationcorrectly describes the Fermi-energy dependence of thesusceptibility in the low-energy region. B. Line-node semimetal
For the line node semimetal Eq. (10), the susceptibilitycan be estimated in an approximate manner as follows.As argued in Sec. II C, the energy band are degener-
FIG. 7. (a) Energy spectrum of the line node HamiltonianEq. (10) with bσ z added ( b ′ = 1 , b = 0 . yz plane. (b) Band structure of thecorresponding lattice model [Eq. (17) with b ′ = 1, b = 0 . u = r = 1] with the slab geometry of the width 100 sites in x direction. The energy bands are plotted against p y , where p z is fixed to 0. ate along the circle node on p y p z -plane. If we regard p z as a parameter and view the Hamiltonian as a two-dimensional system on xy -plane, the energy spectrumon the p x p y -space contains a pair of Weyl nodes when − b ′ /v < p z < b ′ /v . At p z = 0, in particular, the Hamil-tonian can be separated into independent 2 × p z moves away from 0 as shown in Fig. 2(b), we canroughly approximate the whole system by a simplifiedmodel in which the energy spectrum at p z = 0 just per-sists throughout the region of − b ′ /v < p z < b ′ /v .Then the magnetic susceptibility is immediately calcu-lated by integrating the contribution of p z = 0 through − b ′ /v < p z < b ′ /v . Specifically, we substitute Eqs. (34), (35) and (36) with ε sµp z ( x n ) = s p x n + b ,τ sµp z = − s, (46)and X sµp z = X s = ± X µ = ± Z − b ′ /vb ′ /v dp z π ~ , (47)where µ labels the two identical Weyl nodes at fixed p z .As a result, we obtain∆ χ s = e vπ ~ c (cid:16) gµ B c ~ ev (cid:17) b ′ ( −| ε F | )∆ χ so = e vπ ~ c gµ B c ~ ev b ′ θ ( ε F )∆ χ o = − e v π ~ c b ′ δ ( ε F ) . (48)Now the orbital term ∆ χ o has a negative delta-functionas in graphene, and this is natural because the systemis a sum of 2D gapless systems over a finite range of p z .The singularity of the diamagnetism is therefore muchstronger than in the log peak in the point-node case.∆ χ so becomes a step function which has a discontinuityat zero energy. As argued in the previous section, ∆ χ so is formally given by the number of states in the zero-thlevel below the Fermi energy, and now it jumps at ε F = 0because all the zero energy level lies there. The spinmagnetism ∆ χ s linearly decreases as the Fermi energygoes away from zero.The line-node is generally not robust unlike the point-node in the Weyl semimetal, and it can be gapped out bya perturbation. But still, the delta-function diamagneticpeak never just disappears but widens to the gap regionkeeping the total area [i.e. integral of χ ( ε F )], as shownfor the gapped 2D case. We also calculate the susceptibility for the correspond-ing lattice model, Eq. (17). The results are shown asred dashed curves in Fig. 2, where we take u = r = 1and b ′ = 3. As seen from the spectrum in Fig. 2, thecontinuum model is valid only in a small energy range | ε F | < .
5, and the susceptibility inside this region quali-tatively agrees with the approximate expression Eq. (48).Specifically, we see a sharp diamagnetic peak in χ o as wellas the abrupt rise in χ so , while the features are smearedby the finite temperature k B T = 0 .
1. We also see a smallcusp in χ s , which corresponds to the linear decrease inEq. (48).The step-like increase of χ so at zero Fermi energy canbe again understood in terms of the topological sur-face state. Similarly to the point-node case, the finite∆ χ so indicates that the orbital magnetization M o =∆ χ so b/ ( gµ B ) is induced by the spin Zeeman field b . Thestep-like jump in ∆ χ so in Eq. (48) then correspondsto magnetization accumulated near the zero energy, ofwhich amount is M o = eb ′ bπ ~ vc . (49)0If we add the spin field bσ z to the Hamiltonian Eq.(10), the line node is gapped leaving two point nodes at( p y , p z ) = (0 , ± b ′ /v ) as illustrated in Fig. 7(a). In a sys-tem bound by the surfaces perpendicular to x -direction,we have the left and right surface bands to connect thepoint nodes, which are indicated by two green ovals inFig. 7(a). The actual band structure at p z = 0 in the cor-responding lattice model with the slab geometry is shownFig. 7(b), where we see that a pair of surface localizedbands. The surface bands disperse in p y direction and thevelocity is given by v surf = ( b/b ′ ) v . Each surface statecontributes to the electric current ± ev surf counter-flowingthe opposite surfaces. If we adopt the simplification usedabove replacing the spectrum of − b ′ /v < p z < b ′ /v with p z = 0, the area density of the surface electronsis n surf = (2 b ′ /v ) / (2 π ~ ) by considering the momentumregion − b ′ /v < p y , p z < b ′ /v . Finally, the total sur-face current density is j = ev surf n surf = ebb ′ / ( π ~ v ),and the total magnetization given by the surface currentthen becomes M = j/c = eb ′ b/ ( π ~ vc ), which exactlycoincides with Eq. (49). A similar orbital magnetizationinduced by the spin Zeeman field was also argued in thesurface Dirac cone in the topological insulator gapped bythe magnetic perturbation. IV. CONCLUSION
We studied the response to the external magnetic fieldin various three-dimensional nodal systems. We showedthat the magnetic susceptibility consists of the orbitalterm, the spin term and the spin-orbit cross term, andthey have different magnitudes and different character-istic properties. The orbital term is diamagnetic nearthe band touching point, and it exhibits a log divergenceat the point node while a delta-function divergence atthe line node. The spin-orbit cross term is caused bythe spin-orbit interaction, and in the nodal semimetalsit is closely related to the chiral surface current carriedby the topological surface states. In the point-node case,it is a linear function in the Fermi energy while in linenode case, it discontinuously jumps right at the line-node.The spin susceptibility includes the Pauli and Van-Vleckparamagnetism, and in the point node case it is foundto be completely Fermi-energy independent in the linearband regime near the degenerate points. For the calcu-lation, we consider both continuum linear model and thecorresponding lattice model, to show that the continuummodel correctly describes the Fermi-energy dependenceof the low-energy susceptibility except for the constantoffset term.The authors thank Kentaro Nomura for helpful dis-cussions. This work was supported by Grants-in- Aid for Scientific research (Grants No. 25107005).
Appendix A: Susceptibility in general point-nodesemimetals
For the Hamiltonian Eq. (1) with general m and b , thesusceptibility is explicitly calculated using Eqs. (34), (35)and (36) as,∆ χ s = e vπ ~ c (cid:16) gµ B c ~ ev (cid:17) X s,µ θ [ s ( ε F − µb ) − | m | ] × (cid:20) P ( t µ ) + s (cid:18) − ε F + µb (cid:19) t µ (cid:21) , (A1)∆ χ so = e vπ ~ c gµ B c ~ ev X s,µ s θ [ s ( ε F − µb ) − | m | ] t µ , (A2)∆ χ o = e v π ~ c h Q − ( | b | − | ε F | ) θ ( | b | − | m | − | ε F | )+ X µ Q µ ( | ε F | − µ | b | ) θ ( | ε F | − | m | − µ | b | ) i , (A3)where s = ± , µ = ± , and t µ = p ( ε F − µb ) − m ,P ( t ) = 12 t p t + m + m log t + √ t + m | m | ! ,Q µ ( t ) = arccosh tm − µ | b |√ m − b arccos m + µ | b | tm ( t + µ | b | )( | m | > | b | ) , arccosh tm − | b |√ b − m arccosh (cid:12)(cid:12)(cid:12)(cid:12) m + µ | b | tm ( t + µ | b | ) (cid:12)(cid:12)(cid:12)(cid:12) ( | m | < | b | ) . (A4)The magnetization Eq. (32) can also be calculated as M s = e π ~ vc gµ B c ~ ev X s,µ θ [ s ( ε F − µb ) − | m | ] × n ( − µε F + 2 b ) P ( t µ ) + s [ − ε F b + µ ( b + m )] t µ + 13 sµt µ o , (A5) M o = e π ~ vc X s,µ θ [ s ( ε F − µb ) − | m | ] × µ [ P ( t µ ) − s ( ε F − µb ) t µ ] . (A6) S. Murakami, New J. Phys. , 356 (2007). S. M. Young, S. Zaheer, J. C. Teo, C. L. Kane, E. J. Mele,and A. M. Rappe, Phys. Rev. Lett. , 140405 (2012). Z. Wang, Y. Sun, X.-Q. Chen, C. Franchini, G. Xu,H. Weng, X. Dai, and Z. Fang, Phys. Rev. B , 195320(2012). A. Burkov and L. Balents, Phys. Rev. Lett. , 127205(2011). A. Burkov, M. Hook, and L. Balents, Phys. Rev. B ,235126 (2011). X. Wan, A. M. Turner, A. Vishwanath, and S. Y.Savrasov, Phys. Rev. B , 205101 (2011). P. Hosur and X. Qi, Comptes Rendus Physique , 857(2013). Z. K. Liu, B. Zhou2, Y. Zhang, Z. J. Wang,H. M. Weng, D. Prabhakaran, S.-K. Mo, Z. X. Shen,Z. Fang, X. Dai, Z. Hussain, and Y. L. Chen,Science (2014), 10.1126/science.1245085. 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,Nature Materials , 677 (2014). M. Brahlek, N. Bansal, N. Koirala, S.-Y. Xu, M. Neupane,C. Liu, M. Z. Hasan, and S. Oh, Phys. Rev. Lett. ,186403 (2012). S. Borisenko, Q. Gibson, D. Evtushinsky, V. Zabolotnyy,B. B¨uchner, and R. J. Cava, Phys. Rev. Lett. , 027603(2014). S.-Y. Xu, C. Liu, S. K. Kushwaha, R. Sankar, J. W. Krizan,I. Belopolski, M. Neupane, G. Bian, N. Alidoust, T.-R.Chang, et al. , Science , 294 (2015). H. Weng, C. Fang, Z. Fang, B. A. Bernevig, and X. Dai,Phys. Rev. X , 011029 (2015). B. Lv, H. Weng, B. Fu, X. Wang, H. Miao, J. Ma,P. Richard, X. Huang, L. Zhao, G. Chen, et al. , Phys.Rev. X , 031013 (2015). S.-M. Huang, S.-Y. Xu, I. Belopolski, C.-C. Lee, G. Chang,B. Wang, N. Alidoust, G. Bian, M. Neupane, C. Zhang, et al. , Nature communications (2015). S.-Y. Xu, I. Belopolski, N. Alidoust, M. Neupane, G. Bian,C. Zhang, R. Sankar, G. Chang, Z. Yuan, C.-C. Lee, et al. ,Science , 613 (2015). M. Phillips and V. Aji, Physical Review B , 115111(2014). Y. Kim, B. J. Wieder, C. Kane, M. Rappe, et al. , arXivpreprint arXiv:1504.03807 (2015). R. Yu, H. Weng, Z. Fang, X. Dai, and X. Hu, arXivpreprint arXiv:1504.04577 (2015). H. Weng, Y. Liang, Q. Xu, R. Yu, Z. Fang, X. Dai, andY. Kawazoe, Phys. Rev. B , 045108 (2015). S. T. Ramamurthy and T. L. Hughes, arXiv preprintarXiv:1508.01205 (2015). V. Aji, Physical Review B , 241101 (2012). D. T. Son and N. Yamamoto, Physical review letters ,181602 (2012). A. Zyuzin and A. Burkov, Phys. Rev. B , 115133 (2012). M. M. Vazifeh and M. Franz,Phys. Rev. Lett. , 027201 (2013). P. Goswami and S. Tewari, Phys. Rev. B , 245107(2013). D. Son and B. Spivak, Physical Review B , 104412(2013). C.-X. Liu, P. Ye, and X.-L. Qi, Phys. Rev. B , 235306(2013). A. Burkov, Phys. Rev. Lett. , 247203 (2014). J. W. McClure, Phys. Rev. , 666 (1956). S. Safran and F. DiSalvo, Phys. Rev. B , 4889 (1979). M. Koshino and T. Ando, Phys. Rev. B , 235333 (2007). M. Koshino and T. Ando, Physical Review B , 085425(2007). M. Koshino and T. Ando, Phys. Rev. B , 195431 (2010). T. Ito and K. Nomura, Unpublished. Y. Ominato and M. Koshino, Phys. Rev. B , 115433(2013). K. Nomura and N. Nagaosa, Physical Review B , 161401(2010). Y. Tserkovnyak, D. Pesin, and D. Loss, Phys. Rev. B91