Orbital Magnetism of Bloch Electrons I. General Formula

We derive an exact formula of orbital susceptibility expressed in terms of Bloch wave functions, starting from the exact one-line formula by Fukuyama in terms of Green's functions. The obtained formula contains four contributions: (1) Landau-Peierls susceptibility, (2) interband contribution, (3) Fermi surface contribution, and (4) contribution from occupied states. Except for the Landau-Peierls susceptibility, the other three contributions involve the crystal-momentum derivatives of Bloch wave functions. Physical meaning of each term is clarified. The present formula is simplified compared with those obtained previously by Hebborn et al. Based on the formula, it is seen first of all that diamagnetism from core electrons and Van Vleck susceptibility are the only contributions in the atomic limit. The band effects are then studied in terms of linear combination of atomic orbital treating overlap integrals between atomic orbitals as a perturbation and the itinerant feature of Bloch electrons in solids are clarified systematically for the first time.


Introduction
The effect of a magnetic field on electrons in crystals is fascinating and one of the basic problems of solid state physics. Although the fundamental principles are simple, our understanding of this problem has been far from complete due to the complicated matrix elements between different Bloch bands, called as the interband effects of a magnetic field. 1 One of the typical problems of this interband effect is the orbital magnetism.
After the pioneering work of orbital magnetism by Landau for free electrons, 2 the effect of periodic potential was considered by Peierls 3 who obtained the Landau-Peierls formula for the orbital susceptibility, (1.1) where f (ε) is the Fermi distribution function and ε ℓ (k) is the Bloch band energy. This Landau-Peierls formula is obtained by considering the effect of a magnetic field by a phase factor exp iē hc (the so-called Peierls phase) in the hopping integral of the single-band tight-binding model. 3 Here A is a vector potential and e < 0 is the electron charge. Attaching the Peierls phase to the hopping integral corresponds to the modification of the energy dispersion, ε ℓ (k) → ε ℓ (k − eA/ch), in the presence of a magnetic field. Apparently χ LP does not include the deformation of the wave function resulting from the interband matrix elements of * E-mail: ogata@phys.s.u-tokyo.ac.jp the magnetic field. Therefore it is believed that χ LP includes only the intraband effects. (It is not so simple as we show later in the present paper.) Furthermore, χ LP vanishes for insulators since it is proportional to ∂ f (ε ℓ )/∂ ε. On the other hand, in bismuth and its alloys, it has been known experimentally that the diamagnetism takes its maximum when the chemical potential is located in the band gap, [4][5][6][7] i.e., in the insulating state. Apparently the Landau-Peierls formula fails to explain the large diamagnetism in bismuth and its alloys, which had been a mystery for a long time.
Stimulated by this experimental fact, there were various efforts to clarify the interband effects of magnetic field on orbital susceptibility. [8][9][10][11][12][13][14][15][16][17][18][19] For example, several theoretical studies showed that the difference between the total susceptibility and χ LP is in the same order of the difference between χ LP and Landau susceptibility, 9,18,19 even in the nearly-free electron cases. The large diamagnetism of bismuth was finally understood by Fukuyama and Kubo 20 who calculated the magnetic susceptibility based on the Wigner representation. It was clarified that the interband effect of a magnetic field and the strong spin-orbit interaction are essential.
After these theoretical efforts, one of the present authors 21 (hereafter referred as I) discovered an exact but very simple formula of orbital susceptibility as where G is the thermal Green's function G (k, ε n ) in a matrix form whose (i j) component is the matrix element between the i-and j-th band. ε n is Matsubara frequency and γ µ represents the current operator in the µ-direction divided by e/h. The spin multiplicity of 2 has been taken into account and Tr is to take trace over the band indices. Originally this formula is derived based on the Luttinger-Kohn representation. 22 However, as discussed in I, this formula is valid in the usual Bloch representation because the two representations are related by a unitary transformation and the trace is invariant under the unitary transformation. This exact one-line formula (1.3) has been applied to practical models such as Weyl equation realized in graphene and an organic conductor α-(BEDT-TTF) 2 I 3 , [23][24][25] and Dirac equation in bismuth, [26][27][28][29] expressed in the Luttinger-Kohn-type Hamiltonians. For the Bloch representation, the exact formula written in terms of Bloch wave functions had been derived by Hebborn et al. [11][12][13] (especially in Ref. 13 which will be called as HLSS in the following) before the exact one-line formula (1.3) was derived. It was proved in I with the help of the formulation by Ichimaru 17 that (1.3) is equivalent to the results by HLSS. However HLSS's result is very complicated for the practical use. Because of this difficulty, orbital susceptibility for Bloch electrons in general has not been explored in detail. In particular, quantitative estimation of various contributions and their physical meaning have not been clarified.
In the single-band tight-binding model, there is a fundamental problem. When one restricts the band indices of the Green's functions in (1.3) to a single band, one obtains a susceptibility defined as χ 1 which is different from χ LP . 21,[30][31][32] On the other hand, in the two-dimensional honeycomb lattice (or graphene) [30][31][32][33][34][35] which is a typical two-band tightbinding model, it was shown that the orbital susceptibility based on the Peierls phase (eq. (1.2)) is not equal to either χ 1 or χ LP . [30][31][32] From these results, it was claimed that the formula (1.3) cannot be applied to the tight-binding models on one hand, 32 and that there are some "correction terms" to the exact formula (1.3) on the other hand, 30,31 both of which are of course unjustified. As shown in the present paper, these confusions come from the misusage of the exact formula (1.3).
In this paper, starting from the exact one-line formula (1.3) and rewriting it in terms of Bloch wave functions, we derive a new and exact formula of the orbital susceptibility in a different way from those of HLSS. As shown later explicitly, the new formula is equivalent to the previous results. [11][12][13] However, it is simpler than the previous results and contains only four contributions: (1) Landau-Peierls susceptibility, χ LP , (2) interband contribution, χ inter , (3) Fermi surface contribution, χ FS , and (4) contribution from occupied states, χ occ . Except for χ LP , the other three contributions involve the crystalmomentum derivatives of Bloch wave functions. The physical meaning of each term is discussed. In the atomic limit, χ inter is equal to the Van Vleck susceptibility and χ occ is equal to the atomic diamagnetism (or contributions from core-level electrons). χ FS is a newly found contribution proportional to f ′ (ε ℓ ). Then, we apply the present formula to the model of linear combination of atomic orbitals. We will show that the or-bital susceptibility can be calculated systematically by studying the effects of overlap integrals between atomic orbitals as a perturbation from the atomic limit. In this method, itinerant features of Bloch electrons in solids are clarified for the first time. In most of researches, the atomic diamagnetism and Van Vleck contributions are treated separately from χ LP . However, the present exact formula contains all the contributions on the same basis. Furthermore we find that χ occ contains the contributions not only from the core-level electrons (known as atomic diamagnetism), but also from the occupied states in the partially-filled band. This contribution has not been recognized before.
As mentioned above, when we restrict the band indices of the Green's functions in (1.3) to a single band, we do not obtain χ LP . In this paper, we show that the total of several contributions in (1.3) gives χ LP . In these contributions, the fsum rule which involves the summation over the other bands plays important roles. This means that the band indices of the Green's functions in (1.3) should not be restricted to a single band when one consider a single-band tight-binding model, While preparing this paper, we notice that Gao et al 36 studied orbital magnetism in terms of Berry phase using the wavepacket approximation. Their main interest is in the case with broken time-reversal symmetry 37-40 which is not the subject of the present paper. However we can compare our results with theirs in the case where the time-reversal symmetry is not broken. We find that their results are almost equivalent with ours except for a term which has a different prefactor, possibly due to the wave-packet approximation they used.
In section 2 we derive a new formula for orbital susceptibility in Bloch representation. Our main results are summarized in eqs. (2.28)-(2.32) where four contributions are identified. In section 3 we apply the obtained formula to the model of linear combination of atomic orbitals. Section 4 is devoted to discussions and future problems.

Bloch wave functions and current operator
In order to explore the implications of eq. (1.3) in the Bloch representation, some essential ingredients are introduced. Thermal Green's function for the ℓ-th band in (1.3) is simply given by where ε n is Matsubara frequency and ε ℓ (k) is the Bloch band energy. ℓ denotes the band index and the wave vector k is within the first Brillouin zone. In order to obtain the explicit form of the current operator γ µ in (1.3), it is necessary to have information of Bloch wave functions in a periodic potential V (r). From the Bloch's theorem, the eigenfunctions of the Hamiltonian are given by e ik·r u ℓk (r), (2.2) where u ℓk (r) is a periodic function with the same period as V (r) and satisfies the equation Here k · ∇ indicates the inner product between k and ∇. For simplicity, we assume a centrosymmetric potential, V (−r) = V (r). In this case we can choose u † ℓk (r) = u ℓk (−r) 13 where u † ℓk (r) is the complex conjugate of u ℓk (r). For γ µ , we calculate the current operator Substituting the expansion ψ α (r) = ∑ ℓ,kĉ ℓkα e ik·r u ℓk (r), (2.6) into eq. (2.5) and making the Fourier transform, we obtain where the definition of H k in (2.4) has been used. Since γ µ in (1.3) is defined as a current operator divided by e/h, the matrix element of γ µ is given by (see Appendix A) with p ℓℓ ′ µ being the off-diagonal matrix elements 21, 41 Although the integral in (2.9) is sometimes called (interband) "Berry connection", this kind of terms has been familiar for a long time in the literatures. 12,41,42 [Note that the intraband "Berry connection" vanishes in the present Hamiltonian with V (r) = V (−r), as shown in (A·8).] Here we have used the Fourier integral theorem 22 for functions with the lattice periodicity. Originally, the range of the real-space integral on the right-hand side of eqs. (2.7) or (2.9) is within a unit cell, 22 i.e., V Ω Ω · · · d r, where V and Ω are the volumes of the whole system and of the unit cell, respectively. However, the range of integral can be extended to the whole system size V by using the periodicity of u ℓk (r), i.e., V Ω Ω · · · d r = V · · · d r. In the following, the real-space inte-grals are defined in this way.

New formula for orbital susceptibility
Using the above matrix elements for γ µ and thermal Green's functions, we calculate the formula (1.3) in the Bloch representation. Due to the existence of two terms in each γ µ ℓℓ ′ , there appear sixteen terms. Classifying by the band indices of four Green's functions in (1.3), we obtain (2.14) (2.17) where +(x ↔ y) means a term obtained by replacing (x, y) to (y, x). Schematic representation of the contributions to χ 1 -χ 7 are shown in Fig. 1. The summation with prime ∑ ′ means that all the band indices (ℓ, ℓ ′ , ℓ ′′ or ℓ, ℓ ′ , ℓ ′′ , ℓ ′′′ ) are different with each other. In the following we write ε ℓ for ε ℓ (k) as far as it is not confusing. The first contribution χ 1 will be purely intraband since only the intraband matrix elements of γ µ 's are involved. After taking the summation over Matsubara frequency n and making Height of these lines represents the energy level ε ℓ . The array of blue circles connecting the two lines represents the off-diagonal matrix elements of γ µ , i.e., p ℓℓ ′ ,µ . The symbol ⊗ between the two solid lines represents the diagonal component of γ µ , i.e., ∂ ε ℓ /∂ k µ . The right-hand of each diagram is connected to its left-hand because of the trace in (1.3). The red squares indicate the part of the diagrams which can be expressed by the f -sum rule in eq. (2.21).
integrations by parts, we obtain 21 where f (ε ℓ ) is the Fermi distribution function. This χ 1 is similar to the Landau-Peierls susceptibility χ LP in (1.1), but there are two differences. The numerical prefactor of the second term of χ 1 is different from χ LP and the last term of χ 1 does not appear in χ LP . We will show shortly that χ LP is obtained by adding some other contributions from χ 2 , χ 5 and χ 6 . As discussed in Section 1, this means that one should not pick up only χ 1 in discussing the orbital susceptibility in a singleband model. Next let us consider χ 2 . The summation over n in χ 2 gives where the j-th term in χ n is denoted as χ n; j . For χ 2;1 , the summation over ℓ ′ can be carried out and we obtain where we have used the f -sum rule 8, 21 with µ = x, ν = y and the integration by parts. This ℓ ′summation is schematically shown in Fig. 1(b), in which the red square indicates the part of the diagram representing the f -sum rule.
The f -sum rule in eq. (2.21) results from the completeness property of u ℓ ′ k [see eq. (A·13) in Appendix A]. (Various formulas used in the present paper are listed in Appendix A.) One may call the left-hand side of (2.21) as "interband" since it contains the off-diagonal matrix elements of the current operator p ℓℓ ′ ,µ . On the other hand, the right-hand side of (2.21) is expressed by a single-band property, ε ℓ , and as a result, χ 2;1 looks like an "intraband" contribution. This indicates that the naive classification of "intraband" and "interband" does not apply.
For χ 2;2 , the summation over ℓ ′ can also be carried out, and we obtain where we have used The last term of χ 2 is given by (2.24) Here the ℓ ′ -summation can not be carried out due to the presence of the denominator, 1/(ε ℓ − ε ℓ ′ ). This denominator is the same as what appears in the second-order perturbation of the interband process. Features similar to χ 2 are present in χ 5 and χ 6 , as seen in Fig.1, where the excluded terms of ℓ ′ = ℓ ′′ in χ 6 are supplemented by χ 5 , leading to the independent summations over ℓ ′ and ℓ ′′ . As a result, using the f -sum rule, we obtain Now, we can see that sum of χ 1 , χ 2;1 and χ 5;1 + χ 6;1 becomes It is seen that the first two terms give χ LP , while the last term can be combined with other contributions after the transformation which is obtained by putting µντ as xyy in eq. (A·18) in Appendix A.
Other terms in χ 3 -χ 7 are calculated similarly, whose details are shown in Appendix B. We obtain the total susceptibility χ as follows, which is exact as eq. (1.3). (2.32) Schematic representation of these contributions is shown in Fig. 2. It is worth noting that the separation of χ FS and χ occ is not unique. For example, part of χ occ can be rewritten into a form with f ′ (ε ℓ ) by integration by parts. In the above expression, we have chosen a separation which is simple.

Interpretation of each term
The first contribution, χ LP , is Landau-Peierls formula. As shown above, χ LP comes from the contributions of χ 1 , χ 2;1 and χ 5;1 + χ 6;1 . As we have seen, the f -sum rule plays important roles in χ LP .
The second contribution, χ inter , is a purely interband contribution. Among the four contributions, only this term involves two bands of ε ℓ and ε ℓ ′ . In the next section, we show that χ inter is equal to the Van Vleck susceptibility when we consider the atomic limit.
Note that the energy denominator of χ inter is the same as the case of the f -sum rule (2.21). However, there is a clear difference. Let us consider a single-band case in which the energy differences between the bands are quite large. In this case, we have ∆E = min|ε ℓ (k) − ε ℓ ′ (k)| which satisfies ∆E >> [band width of ℓ-th band]. Then, the absolute value of χ inter is less (2.33) Then the ℓ ′ summation can be carried out using the completeness property of u ℓ ′ k . For example, a typical term can be written as which is arbitrarily small when ∆E is large. This is in a sharp contrast with the f -sum rule (2.21) in which the right-hand side is finite when ∆E is large. We call the third contribution in (2.28) as χ FS (FS stands for "Fermi surface"), since it is proportional to f ′ (ε ℓ ). This is a newly found contribution, but its physical meaning is not clear at present. However, the factor in the integral in χ FS is common with χ inter , indicating some close relationship between χ FS and χ inter . The fourth contribution, χ occ , has contributions from the occupied states ("occ" stands for occupied states). As shown in the next section, χ occ is equal to the atomic diamagnetism in the atomic limit. Furthermore, we show that χ occ contains the contributions not only from the core-level electrons, but also from the occupied states in the partially-filled band. (See the right-hand side of Fig. 2.) This contribution has not been recognized before.

Comparison with the result by HLSS
For the orbital susceptibility in (2.28), we have only four contributions which are simpler than those obtained previously by HLSS, i.e., eqs. (4.3)-(4.6) in Ref.. 13 In Appendix C, we prove the equivalence between the present result and HLSS. Here we summarize the differences between the two.
(1) The result by HLSS is not symmetric with respect to the exchange of x and y. This is because they used Landau gauge, A = (−Hy, 0, 0). On the other hand, the present formula is symmetric with respect to x and y because we have used the gauge-invariant formalism (1.3). In order to prove the equivalence between our result and that by HLSS, we have to symmetrize the HLSS's result. (See details in Appendix C.) (2) Among the four contributions in the present formula, χ LP is determined solely from the energy dispersion ε ℓ (k).

The other three contributions involve the k-derivatives of
Bloch bands (See eq. (C·1).) As shown above in eq. (2.27), this term can be rewritten in terms of u ℓk and has been included in χ FS in our formalism. It is important to use (2.27) in order to simplify the final expression. Note that, in contrast to (2.36), χ LP can not be rewritten in terms of u ℓk 's.
(3) The result by HLSS contains several terms which have a common denominator of 1/(ε ℓ − ε ℓ ′ ). In our formula, these contributions are summed up into a single term as χ inter . The method how we can sum up the several terms in HLSS into a single term is explained in Appendix C.
(4) As explained above, each contribution in the present formula has a rather clear meaning compared with the previous ones. For example, χ inter and χ occ are contributions naturally connected to the Van Vleck susceptibility and atomic diamagnetism, respectively. Note that, in the HLSS's formula, the contribution of the atomic diamagnetism is expressed as the first term of χ However, the numerical prefactor is different from that of the present result. (See the term proportional toh 2 /m in χ occ ). As shown in the next section, this term reproduces the orbital susceptibility from the core-electrons in the atomic limit. χ occ in the present formula gives a correct prefactor, while eq. (2.37) by HLSS does not. As shown in Appendix C, we find that the correct term is obtained in the HLSS's formula when we rewrite the interband contributions into a single term as χ inter .

Comparison with the results by Gao et al.
Let us discuss here the recent work by Gao et al 36 who studied orbital magnetism in terms of Berry phase. They are interested in the case with broken time-reversal symmetry in which spontaneous orbital magnetization appears. [37][38][39][40] In this case, there are several terms involving the Berry curvature denoted as Ω. In the present notation, its z-component is given by However, in our case with a centrosymmetric potential, Ω z vanishes since there is a relation u † ℓk (−r) = u ℓk (r). As a result, we do not have the contributions coming from the Berry curvature. However, we can compare our results with theirs in the case where the time-reversal symmetry is not broken.
Details of calculations are shown in Appendix D. By using the completeness property of u ℓk , we can show that their results are almost equivalent with our results except for the coefficient of χ FS . We think that this difference is due to the wave-packet approximation 39 used in their formalism. Nevertheless, their formula for the orbital susceptibility based on the wave-packet approximation is fairly accurate.

Band effect from atomic limit
The obtained formula in eqs. (2.29)-(2.32) is exact. However, in order to calculate each contribution explicitly, it is necessary to specify the functional form of u ℓk (r). For example, u ℓk (r) can be obtained in general from the first-principle band calculation. In this paper, however, we will study each contribution from the atomic limit to see possible band effects based on the linear combination of atomic orbitals (LCAO). In the atomic limit, it is found that diamagnetic susceptibility from core electrons and Van Vleck susceptibility are the only contributions to χ. Then, χ is estimated by treating the overlap integrals between atomic orbitals as a perturbation. We will show that there appear several contributions to χ in addition to χ LP . In this perturbative method, the itinerant feature of Bloch electrons in solids are clarified systematically. Schematic picture is shown in Fig. 2.

Atomic limit
In order to study the atomic limit in the present formula, it is appropriate to use LCAO. Let us consider a situation in which the periodic potential V (r) is written as where R i represent lattice sites and V 0 (r) is a potential of a single atom. We use atomic orbitals φ n (r) which satisfy −h 2 2m ∇ 2 + V 0 (r) φ n (r) = E n φ n (r).

(3.2)
Using these atomic orbitals, we consider the LCAO wave function which is used as a basis set for u ℓk (r). Here N is the total number of unit cells. It is easily shown that ϕ nk (r) are periodic functions with the same period with V (r).
In the atomic limit, V 0 (r − R i ) and φ n (r − R i ) are confined in a unit cell and there is no overlap between nearest-neighbor V 0 (r − R i ) or nearest-neighbor φ n (r − R i ). In this case, it is easily shown that the LCAO wave function, ϕ ℓk (r), in (3.3) satisfies the equation (2.3) with energy ε ℓ = E ℓ that is independent on k. Therefore, u ℓk is just given by By substituting eq. (3.4) and ε ℓ = E ℓ into eqs. (2.29)-(2.32), we obtain χ LP , χ inter , χ FS and χ occ . Since E ℓ is k-independent, χ LP = χ FS = 0. For χ inter , we obtain Since there is no overlap between atomic orbitals, only the terms with R i = R j survives. Thus χ inter is simplified as where the r-integral is shifted to the center of the atomic orbital, R i . The right-hand side of (3.6) is nothing but the Van Vleck susceptibility which we denote as χ (Van Vleck) . Similarly, by substituting (3.4) and ε ℓ = E ℓ into (2.32), we This is just the atomic diamagnetism coming from the core electrons which we denote as χ (atomic dia.) . Therefore, in the atomic limit, we have χ = χ (Van Vleck) + χ (atomic dia.) . (See the left-hand side of Fig. 2.)

Perturbation with respect to the overlap integrals
Next we consider the case in which there are overlap integrals between the nearest-neighbor atomic orbitals. In this case, using ϕ nk (r) in eq. (3.3), we expand u ℓk as where the Hamiltonian matrix elements are h nm (k) = ϕ * nk (r)H k ϕ mk (r)d r, (3.10) and s nm (k) represents the integral s nm (k) = ϕ * nk (r)ϕ mk (r)d r. (3.11) h nm (k) and s nm (k) can be calculated perturbatively with respect to the overlap integral (3.12) with O being an operator and R j R i . For example, the firstorder term of h nm (k) contains the hopping integral used in the tight-binding model. By substituting eq. (3.8) and ε ℓ (k) into eqs. (2.29)-(2.32), we can show that each of four contributions, χ LP , χ inter , χ FS and χ occ , is calculated perturbatively with respect to the overlap integrals. In contrast to the atomic limit, there are two new features: (1) ε ℓ (k) has band dispersion due to the hopping integrals, and (2) u ℓk (r) has an additional k-dependence through c ℓ,n (k) in eq. (3.8). The latter gives several contributions to the orbital susceptibility originating from the kderivatives of u ℓk (r). One may expect that χ LP is dominant in the first-order perturbation. However, we find that the situation is not so simple even in the single-band case. Each of χ LP , χ inter , χ FS and χ occ depends on the location of the chemical potential as well as on the details of the model. In the forthcoming paper, we will discuss several explicit models such as single-band and two-band tight-binding models.

Discussions
Based on the exact formula, we have shown rigorously that the orbital susceptibility for Bloch electrons can be described in terms of four contributions, χ = χ LP + χ inter + χ FS + χ occ . Except for the Landau-Peierls susceptibility, χ LP , the other three contributions involve the crystal-momentum derivatives of u ℓk 's. These contributions represent the effects of the deformation of the wave function due to the magnetic field. We find that χ occ contains the contributions from the occupied states in the partially-filled band, which has not been recognized before. We applied the present formula to the model of LCAO. In the atomic limit where there are no overlap integrals, χ inter becomes χ (Van Vleck) and χ occ becomes χ (atomic dia.) . These two are the only contributions to χ in the atomic limit. When the overlap integrals are finite, we have discussed that χ can be calculated by treating the overlap integrals as a perturbation. In this method, itinerant features of Bloch electrons in solids can be clarified systematically for the first time.
The present formalism can be used as a starting point for various extensions. Several future problems are as follows: (1) It is very interesting to apply the present formula to the multi-band tight-binding models. A typical example is the honeycomb lattice which is a model for graphene. [30][31][32][33][34][35] In this case, we have A-and B-sublattice in a unit cell, and as a result, we have massless Dirac electrons (or more precisely, Weyl electrons) which is a typical two-band model. The orbital susceptibility has been calculated by several groups 30-32 based on the Peierls phase. In contrast, in the present formula, all the contributions from Bloch bands are included rigorously. Application of the present formula to graphene will be discussed in the forthcoming paper.
(2) In the present Hamiltonian, the spin-orbit interaction is not included. It is also a very interesting problem to study the orbital susceptibility in the presence of spin-orbit interaction. As discussed recently by Gao et al., 36 the orbital susceptibility in the Hamiltonian with broken time-reversal symmetry [37][38][39][40] is another interesting problem. This will be also studied in the forthcoming paper. As suggested 36 there appear several terms which is written in terms of Berry curvatures.
(3) We have confined ourselves in the orbital susceptibility in this paper. The transport coefficients are of course interesting quantities. 43 Hall conductivity in the Weyl equation realized in graphene and an organic conductor α-(BEDT-TTF) 2 I 3 23-25 as well as in bismuth 26,29 has been discussed. The similar method used in this paper can be applied to the Hall conductivity in the Bloch representation.

Acknowledgment
We thank very fruitful discussions with F. Piéchon, I. Proskurin, Y. Fuseya, H. Matsuura, T. Mizoguchi and N. Okuma. This paper is dedicated to Professor Ryogo Kubo (1920Kubo ( -1995, who guided authors into the never-fashionable but very deep and fundamentally important problem of orbital magnetism in solids. This work was supported by a Grantin-Aid for Scientific Research on "Dirac Electrons in Solids" (No. 24244053) and "Multiferroics in Dirac electron materials" (No. 15H02108).

Appendix A: Various formulas for matrix elements of Bloch wave functions
In this appendix, we derive various useful formulas which are quite often used in the derivations. First, we make k derivative of the equation for u ℓk in eq. (2.3) In the following, we write ε ℓ for ε ℓ (k). When we multiply u † ℓ ′ k and make the real-space integral, we obtain Exchange of ℓ ↔ ℓ ′ gives eq. (2.8).
Next, the k µ and k ν derivative of eq. (2.3) gives When we multiply u † ℓ ′ k and make the real-space integral, we obtain (A·4) Similarly, when we multiply ∂ 2 u † ℓk /∂ k τ ∂ k σ and make the real-space integral, then we obtain We can obtain other series of formulas by making k-derivative of the ortho-normal condition, u † ℓ ′ k u ℓk d r = δ ℓℓ ′ . First we obtain Furthermore if the system is centrosymmetric, i.e., if V (−r) = V (r) holds, we can choose u † ℓk (−r) = u ℓk (r). 13 In this case, we can make further convenient formulas. Firstly using the change of variable r → −r, we obtain Combining this relation with (A·6) with ℓ = ℓ ′ , we obtain (This integral is the so-called intraband "Berry connection".) This formula is used quite often. Further k-derivatives give Using the relation u ℓk (−r) = u † ℓk (r), we can also show that p ℓℓ ′ µ in (2.9) is real and Furthermore, from the k ν derivative of eq. (A·2) with ℓ = ℓ ′ , we obtain Note that µ and ν can be exchanged (µ ↔ ν) in (A·11). Finally, using (A·11) and (A·1), we have By using p ℓℓ ′ ,µ in (2.9), we obtain the f -sum rule 8,21 as follows: where we have added the ℓ ′ = ℓ term which is equal to zero due to the factor (ε ℓ − ε ℓ ′ ). We have also used the relations (A·1) and (A·6), and the completeness property, ∑ ℓ ′ u ℓ ′ k (r)u † ℓ ′ k (r ′ ) = δ (r − r ′ ). Similarly we obtain where we have added the ℓ ′′ = ℓ term in (A·14) and ℓ ′ = ℓ term in (A·15) which are equal to zero since u † ℓk (∂ u ℓk /∂ k µ )dr = 0 (eq. (A·8)).
When we multiply u † ℓk to the k µ , k ν , k τ derivative of eq. (2.3) and make the real-space integral, we obtain (A·16) where we have used (A·8). On the other hand, k τ derivative of (A·11) gives Substituting this relation into the right-hand side of (A·16) and using (A·9), we can show 1 2 (A·18)