Orbital Magnetism of Bloch Electrons III. Application to Graphene

The orbital susceptibility for graphene is calculated exactly up to the first order with respect to the overlap integrals between neighboring atomic orbitals. The general and rigorous theory of orbital susceptibility developed in the preceding paper is applied to a model for graphene as a typical two-band model. It is found that there are contributions from interband, Fermi surface, and occupied states in addition to the Landau--Peierls orbital susceptibility. The relative phase between the atomic orbitals on the two sublattices related to the chirality of Dirac cones plays an important role. It is shown that there are some additional contributions to the orbital susceptibility that are not included in the previous calculations using the Peierls phase in the tight-binding model for graphene. The physical origin of this difference is clarified in terms of the corrections to the Peierls phase.


Introduction
Graphene has a simple and interesting electron system that contains massless chiral Dirac electrons (or Weyl electrons) in a two-dimensional honeycomb lattice. Its various electronic properties have been explored extensively. [1][2][3][4] Among them, orbital magnetism is an interesting property since the Dirac electrons have strong interband effects 5) between the upper and lower Dirac cones.
Actually in early research, McClure 6) showed that orbital susceptibility has a delta function-like peak as a function of chemical potential µ where the two Dirac cones come in contact with each other. This is confirmed by using the exact oneline formula (Fukuyama formula). 7,8) However, these calculations assume a linear dispersion ε = ±v|k| with a finite cutoff. In actual graphene, on the other hand, the energy dispersion deviates from the linear dispersion away from the Dirac points in the Brillouin zone. Therefore, it is necessary to take account of the Bloch bands. Although there have been several studies on the magnetic susceptibility for graphene, [9][10][11][12][13][14] we revisit this issue in the present paper by performing a systematic expansion with respect to the overlap integrals between nearest-neighbor atomic orbitals based on an exact formula expressed in terms of Bloch wave functions and the energy dispersion. 15,16) The orbital susceptibility for graphene or a twodimensional honeycomb lattice was studied [9][10][11] using the Fukuyama formula 7) Tr γ x G γ y G γ x G γ y G , (1.1) where G represents the thermal Green's function G (k, ε n ) of the matrix form with respect to the band indices, ε n is the Matsubara frequency, and γ µ (µ = x, y) is the current operator in the µ-direction divided by e/h. The spin multiplicity of 2 has been taken into account and Tr denotes the trace over the band indices. This formula is exact if all the bands are taken into account. However, in the calculations, [9][10][11] the band indices of the Green's functions were restricted to the upper and lower Dirac cones. Since the band indices of the Green's functions in Eq. (1.1) should not be restricted to a few bands, as discussed in detail in Ref. 15, the results should be reexamined.
Recently, by claiming that there are some "correction terms" 12,13) to the exact formula in (1.1) in the case of the two-band model, the orbital susceptibility for graphene was studied. Calculations based on the Peierls phase in the tightbinding model 14) also gave the same susceptibility. More recently, Gao et al. 17) studied a model for gapped graphene based on the wave-packet formalism and obtained consistent results with those obtained from the Peierls phase. 14) Although the above results are consistent with each other, we examine this issue motivated by the following findings. Recently we studied the orbital susceptibility for single-band models based on an exact formalism 15,16) (referred to as I and II in the following). It was found that there are comparable contributions in addition to the susceptibility originating from the Peierls phase. 16) Furthermore, we clarified the corrections of the Peierls phase in the tight-binding model. 18) Therefore, we expect that there are additional contributions also in graphene. For this purpose, we think it is important to calculate the susceptibility by systematic expansion with respect to the overlap integrals, as carried out for single-band models in II. 16) In our preceding paper I, 15) we rewrote the Fukuyama formula in (1.1) in terms of Bloch wave functions and obtained a new and equivalent formula for the orbital susceptibility χ as follows: (1. 2) The suffixes of χ LP , χ inter , χ FS , and χ occ denote Landau-Peierls, interband, Fermi surface, and occupied states, respectively. 15) Note that the formula in (1.2) is exact, as is Eq. (1.1). As shown in II, 16) χ LP is in the first order with respect to overlap integrals. Thus, we calculate each term in (1.2) for graphene up to the same order with χ LP . In contrast to the previous studies, [9][10][11][12][13][14] this is the first exact calculation up to the first order with respect to the overlap integrals starting from the atomic limit. We will show that there are some contributions that were not included before. The physical origin of these additional contributions is discussed in terms of the corrections to the Peierls phase in the tight-binding model. 18) In the single-band model studied in II, χ inter vanishes. 16) In contrast, in graphene, which is a typical two-band model, we show that χ inter contributes to the total susceptibility. This paper is organized as follows. In Section 2, we develop an atomic orbital model for graphene. We develop a usual tight-binding model but the wave functions are explicitly obtained in terms of p π atomic orbitals for carbon atoms. Then we calculate the orbital susceptibility in Section 3 by systematic expansion with respect to the overlap integrals. It is found that the relative phase between the atomic orbitals for the A and B sublattices plays an important role leading to contributions comparable to χ LP . Section 4 is devoted to discussion and summary.

p π orbitals and Hamiltonian
Since there are two carbon atoms in a unit cell (A and B sublattices) as shown in Fig. 1, we have two bands that form massless Dirac electrons, or Weyl electrons. First, we construct the p π (or p z ) band for graphene. [In the following, we do not consider the contributions from the core-level electrons, i.e., the 1s orbital, or the 2s, p x , and p y orbitals forming σ -bonds. We focus on the contributions from the p π orbitals of carbon atoms considering that only the p π band crosses the Fermi energy. ] As in I and II, 15,16) we assume that the periodic potential V (r) is written as where R i represents the positions of carbon atoms forming a honeycomb lattice. Under the periodic potential V (r), wave functions are given by e ik·r u ℓk (r), where u ℓk (r) satisfies For the carbon 2p π orbital, it was discussed that the nucleus charge is screened by core electrons, and as a result, the effective atomic potential is given by where Z eff represents the effective charge, Z eff = 3.25. 19) The p π (or p z ) atomic orbital is then given by where a * B is the renormalized Bohr radius with a B =h 2 /me 2 . As in II, 16) we use the orthogonal wave functions 20) where the j-summation represents the sum over the nearestneighbor (n.n.) sites of R i and s is the overlap integral Note that s is independent of the direction R = R j − R i since the p π orbital is isotropic in the xy-plane. In the following, we calculate the orbital susceptibility up to the first order with respect to "overlap integrals" whose integrand contains the overlap of atomic orbitals, φ * Using these orthogonal wave functions, we consider two linear combinations of atomic orbitals (LCAOs) for the A and B sublattices defined as and ϕ ortho Here N is the total number of sites on each sublattice and R Ai (R Bi ) represents the position of the site in the A (B) sublattice in the i-th unit cell. We assume that u ℓk (r) can be expanded in terms of ϕ ortho Ak (r) and ϕ ortho Bk (r). The mixing of other orbitals is neglected. In this approximation, the matrix elements of H k are given by where n, m = A or B. Considering that the nearest-neighbor sites are A and B sublattices, we obtain 16)  where E pπ is the atomic energy eigenvalue for the p π orbital and (2.13) The derivations of t 0 and c pπ in the hopping integral t were discussed and justified in II. 16) γ k is given by where R are the vectors from a B site to its nearest-neighbor A sites, R = R A j − R Bi . (Here, we have assumed only the nearest-neighbor hopping integrals, but the extension to longer-range hopping integrals is straightforward.) Explicitly, γ k becomes γ k = e −ik y a + e i( where a is the distance between the nearest-neighbor sites, i.e., a = |R|. For the p π -orbital, the integrals can be calculated explicitly as 21) Note that these integrals do not depend on the direction of R. s and t are in the first order with respect to the overlap integrals, which means that they are proportional to e −p . Using a = 1.42A in graphene, we obtainp = 4.36, s = 0.237, and t = 3.55eV. The overlap integrals and various matrix elements are shown in Appendix A.

Energy dispersion and chirality around Dirac points
By diagonalizing the Hamiltonian in (2.12), we obtain normalized eigenfunctions and energy eigenvalues as follows: and is the antibonding (bonding) state between the two sublattices with phase factor e ±iθ k /2 . ε + pπ (k) (ε − pπ (k)) gives the energy dispersion of the upper (lower) Dirac cone. As in II, 16) the constant energy E pπ + C pπpπ is included in the chemical potential in the following, and we write the Fermi distribution function The phase factors e ±iθ k /2 in u ± pπ k (r) are determined in order to satisfy These relations are required in the case of a centrosymmetric potential and have been used in various steps to derive Eq. (1.2) in I. 15) To prove Eq. (2.20), it is necessary to note that each A site at R Ai has its partner of the B site satisfying −R Ai = R B j . Using this relation, we can see that ϕ ortho Ak (−r) = ϕ ortho * Bk (r) holds, which leads to Eq. (2.20). The density of states for the honeycomb lattice was obtained in the context of the phonon density of states. 22) It is given by where K(κ) is the elliptic integral of the first kind with Note that E pπ + C pπpπ is included in µ. For completeness, its derivation and some limiting cases are given in Appendix B. Two Dirac points exist at k = k ± 0 = (± 4π 3 √ 3a , 0), and around these Dirac points, γ k can be expanded as where η is the angle between the k x axis and vector k. Thus, the phase θ k of γ k is related to the chirality for each Dirac point, and its chirality is opposite for the two Dirac points.
In the following calculations, we find that the k-derivatives of θ k play important roles since θ k is attached to the wave function u ± pπ k (r). Furthermore, we find that ∂ θ k /∂ k µ (µ = x, y) is related to an integral as as shown in Appendix A. This integral can be called the interband "Berry connection" between the upper and lower Dirac cones. However, this kind of integral has been familiar for a long time in the literature. [23][24][25]

Orbital susceptibility for Graphene
Using the energy dispersion in the previous section, the Landau-Peierls susceptibility 26,27) from the p π orbital is given by where we have used the abbreviations which is in the first order with respect to the overlap integrals. As a result, χ LP is also in the first order of the overlap integrals. In the following, we calculate χ inter , χ FS , and χ occ up to the same order.
To evaluate χ inter , χ FS , and χ occ , we use with C + = 1, C − = −i and we have also used abbreviations such as θ x = ∂ θ k /∂ k x . In the present model, χ inter is given by 15,16) where the range of the real-space integral · · · d r has been extended to the whole system size by using the periodicity of u ℓk (r), 15) and (x ↔ y) represents terms in which x and y are exchanged. Here, we consider the cases where ε ℓ = ε + pπ (k) and ε ℓ = ε − pπ (k). Because of the Fermi distribution function f (ε ℓ ), these cases give the contributions from the occupied p π orbitals. [As discussed before, we do not consider the contributions from the other occupied bands, i.e., from the 1s, 2s, p x , and p y orbitals.] In contrast, we need to take the summation over ℓ ′ (ℓ ′ ℓ) in Eq. (3.4), because these terms originate from the virtual processes in the second-order perturbation.
In the case of the 1s single-band model discussed in II, 16) the second line of (3.4) vanishes both in the zeroth and first order with respect to the overlap integrals sinceL z φ 1s (r) = 0. However, in the present two-band model, there is a contribution in the zeroth order that involves the derivatives of the phase θ k even ifL z φ p π (r) = 0. This is in sharp contrast to the single-band case. Since there is a zeroth-order term, the infinite summation of ℓ ′ in χ inter should be carried out carefully to obtain the contributions up to the first order of the overlap integrals.
Furthermore, it should be noted that the denominator in Since this denominator is in the first order with respect to the overlap integrals, this case should be treated carefully. As shown in Appendix C, however, the numerator in this case turns out to be proportional to the fourth order with respect to the overlap integrals. Therefore, the perturbation does not break down. In this case, we find that the combination of ∂ H k ∂ k y + ∂ ε ℓ ∂ k y in Eq. (3.4) is important. As shown in Appendix C, the summation over ℓ ′ in (3.4) is carried out analytically, and we obtain Note that all the terms are related to the derivatives of θ k . In the single-band case, 16) the relative phase θ k does not appear and thus χ inter vanishes. Next we calculate χ FS from the p π orbital, which is given by 15,16) (3.6) These integrals can be rewritten with the help of the relation L z φ p π (r) = 0 in a similar way to χ inter . Then we obtain (3.7) Details of the derivation are shown in Appendix C. Finally, χ occ is given by (3.8) [Here again, we do not consider the contributions from the other occupied bands.] In a similar way carried out in our pre-ceding paper, 16) χ occ becomes where Re denotes the real part and the expectation value is defined as Note that there is a difference between Φ pπ (r) and φ pπ (r) as in Eq. (2.7). Using the expectation values in (A·2), we obtain where we have used the definition of γ k in (2.14). We numerically calculate χ LP , χ inter , χ FS , and χ occ as a function of chemical potential µ at T = 0. The results are shown in Fig. 2, in which each contribution is normalized by the Pauli susceptibility χ 0 at the band edge (µ = ±3t), i.e., where the system size is L 2 = 3 √ 3a 2 N/2 with N being the total number of unit cells. Here, we have used the fact that the model is equivalent to free electrons with effective mass m * = 2h 2 /3ta 2 at the bottom of the band.
There are several remarks on the above results.
(1) In contrast to the single-band case discussed in I, 15) there appear several terms involving θ x and θ y , which originate from the two-band nature. Furthermore, χ inter is nonzero, which is in sharp contrast to the single-band case, where χ inter vanishes. This is also owing to the two-band nature.
(2) There are three terms that are in the zeroth order with respect to the overlap integrals: the first term of χ inter in Eq. (3.5) and the first two terms of χ occ in Eq. (3.11). However, the first term of χ inter and the second term of χ occ , both of which diverge as µ → 0, exactly cancel with each other. Therefore, only the first term of χ occ contributes to the total susceptibility in the zeroth order. This term is proportional to the electron  number in the p π band, i.e., where n(µ) represents the total electron number with spin degeneracy when the chemical potential is µ, and it is calculated from the density of states in (2.21) as n(µ)/L 2 = 2 µ −3t D(µ)dµ. χ occ:1 represents the contributions from the occupied states in the partially filled p π -band, which we call intraband atomic diamagnetism. 16) Since χ occ:1 does not have a factor of e −p , it gives comparative contributions as χ LP in a similar way to the single-band case studied in II.
(3) At the band bottom (µ = −3t), only χ LP has a contribution. Its value is just equal to −1/3χ 0 , which is understood as Landau's diamagnetic orbital susceptibility for free electrons. Furthermore, χ LP has a diverging peak at µ = −t, which corresponds to the van Hove singularity.
(4) χ FS is always negative and its small wiggle at µ = −t is owing to a subtle cancellation between the last two terms in (3.7), both of which diverge as µ → 0. χ FS has a sizable contribution in the region of −3t < µ < −t.
The total susceptibility χ = χ LP + χ inter + χ FS + χ occ is shown in Fig. 3 as a function of µ. A δ -function-like peak at µ = 0 6) is not included, which we discuss shortly. It is rather surprising that the total of each contribution becomes a smooth function of µ irrespective of the irregular µ-dependences of χ inter and χ FS near µ = −t in Fig. 2. In  Fig. 3, the present result is compared with that obtained by Gòmez-Santos and Stauber, 13) denoted as χ (GS) . Apparently, there is a sizable difference from χ (GS) , which is discussed in the next section.
In the present results, there is asymmetry with respect to the sign change of µ. This is because the contribution χ occ:1 in Eq. (3.13) is a monotonically decreasing function of µ. When the p π band is fully filled (i.e., µ > 3t), only χ occ:1 gives the contribution to the orbital susceptibility. Therefore, the total susceptibility becomes 14) with N e being the total electron number of the p π band. This is nothing but the atomic diamagnetism from the p π electrons.

Discussion and Summary
First, in order to see the relationship between the present result and the previous ones, we rewrite the total susceptibility in a different way. Using integration by parts for χ FS in (3.7), we find that the total susceptibility can be rewritten as the following simple form: Note that only χ occ:1 is in the zeroth order with respect to the overlap integrals, while the other three contributions are in the first order. χ occ:1 originates from χ occ , but χ 1 and χ 2 are combinations of χ inter , χ FS , and χ occ . Finally, when we use the interesting relation (see Appendix D), χ 2 can be rewritten as where we have used the relation in (A·5). Expectation values for an operator O in terms of φ pπ (r) are defined as (4.7) Using the expectation values calculated in Appendix A, each term in Eq. (4.6) becomes 3a 2 a * 2 B = 0.0395a 4 , (4.8) whenp = a/2a * B = 4.36. Each contribution, χ occ:1 , χ LP , χ 1 , and χ 2 , is plotted in Fig. 4 together with the total susceptibility. The flat part of the total susceptibility near µ = 0 originates from the cancellation between the µ-dependences of χ LP and χ 1 . As discussed before, only χ occ:1 has asymmetry with respect to ±µ.
Here we comment on the δ -function-like peak at µ = 0, 6) which is not included in Figs. 3 and 4. As shown by Fukuyama,8) if the energy dispersion is completely linear around a Dirac point, i.e., ε = ±v|k|, we expect as a contribution from the Dirac cones at T = 0. Here Γ is phenomenologically introduced damping and the presence of two Dirac points has been taken into account. This result at T = 0 with finite damping 8) corresponds to the δ -function-like peak obtained by McClure 6) in clean systems at finite temperatures. This peak will originate from the singular behavior of the kderivatives of θ k at the Dirac points. In the present model, we have v = 3ta/2 from Eq. (2.23) and thus χ (Dirac) is proportional to t 2 , which means that χ (Dirac) will appear in the second order with respect to the overlap integrals. Although χ (Dirac) is in the second order, it should be included in the total susceptibility because of its singular behavior. Let us compare the present result with previous ones. Gòmez-Santos and Stauber 13) used the following formula for orbital susceptibility: Tr γ x Gγ y Gγ x Gγ y G where G is now a 2 × 2 matrix, G = (iε n − H k ) −1 , andγ µ = ∂ H k /∂ k µ , etc., with H k being the Hamiltonian in a 2 × 2 matrix form. (The spin degeneracy has been included.) The second term in (4.10) is the "correction term" 13) added to the exact formula of (1.1). Actually, before Gòmez-Santos and Stauber, Koshino and Ando 12) used another formula, although they did not calculate χ (KA) explicitly. The last three terms are their "correction terms". We can show that χ (GS) and χ (KA) are equivalent using the relation ∂ G k /∂ k µ = Gγ µ G and integration by parts.
In the present notation, the 2 × 2 Hamiltonian is given by with σ x,y being the Pauli matrix. After some algebra, we can show that χ (GS) is given in the present notation as ± f (±ε k ) (ε x θ y θ xy + ε y θ x θ xy + ε xy θ x θ y ) .

(4.13)
When we perform integration by parts in the second summation, we can see that χ (GS) is exactly equal to χ LP + χ 1 in Eq. (4.1). Therefore, the difference between the present result and χ (GS) is simply χ occ:1 and χ 2 . The origin of this difference will be the following. In χ (GS) , (i) the effect of deformation of the wave function u ℓk is not taken into account completely and (ii) the occupied electron contributions in χ occ are not included. (iii) As shown in I, 15) the important interband contributions are taken into account through the f -sum rule. However, "the correction terms" introduced in the previous studies do not sufficiently take account of these contributions.
In χ (GS) and χ (KA) , the effect of a magnetic field was introduced in the energy dispersion of the Bloch bands. The correct procedure should be to introduce the effect of a magnetic field into the Bloch equation in (2.2). As a result, the exact Fukuyama formula in Eq. (1.1) is obtained, which contains all the contributions including the effects of deformation of the wave functions.
As explained in the introduction, calculations based on the Peierls phase 14) in the tight-binding model give the same susceptibility as χ (GS) . The equivalence of these two methods is shown in Appendix E. This result means that the calculations based on the Peierls phase do not contain χ occ:1 and χ 2 . Therefore, we expect that the effects of a magnetic field other than the Peierls phase will play important roles. Recently, we clarified corrections to the Peierls phase argument in the tightbinding model and studied the cases of a benzene molecule and a square lattice. 18) A similar argument can be applied to graphene. By extending Pople's formulation on the effects of a magnetic field to the tight-binding model, 28) we showed that the hopping integral t has the correction ∆t = t 2h 2 with h = eHa 2 /2ch. 18) (The expression for t 2 is shown shortly.) This correction originates from the deformation of the wave function by the magnetic field. Since a change in the hopping integral causes a change in the kinetic energy, it contributes to the susceptibility as This formula is equivalent to χ 2 in Eq. (4.5), and b in (4.5) corresponds to a 4 t 2 /t. In the present notation, t 2 is given by 18) When we numerically evaluate the expectation values, we find that a 4 t 2 /t is slightly different from b. This is probably owing to the difference in the evaluation of integrals involving V 0 (r − R), and we think that this difference is not important. Therefore, ∆χ qualitatively corresponds to χ 2 . To confirm this correspondence, we estimate the second term in (4.15) as follows: where we have used the definition of t in (2.13). Using the expectation value in (A·2) and R 2 x + R 2 y = a 2 , (4.16) becomes In this estimation, t 2 becomes (4.18) Substituting this expression into (4.14), we recover the results in (4.5) and (4.6).
Finally, let us discuss the relationship of the results of this work to the recent work by Gao et al. 17) [Note that their formula in the time-reversal symmetric case is almost equivalent to the present formula in (1.2) except for the numerical factor of χ FS . 15) ] In Ref. 17, the orbital susceptibility for gapped graphene was calculated. It was claimed that the contribution of energy polarization (which corresponds to our χ FS /2) and that of the Van Vleck susceptibility (which corresponds to χ inter ) vanish in the case of gapped graphene. Furthermore, Langevin-type and geometrical susceptibilities (which correspond to χ occ ) give a symmetric contribution with respect to the sign change of µ. These results are different from the present case of (gapless) graphene. Thus, it is an interesting future problem to examine the model of gapped graphene, in which the potential V (r) is noncentrosymmetric. In such cases, orbital magnetization 29) and the connection to the socalled "Berry" curvature can be discussed.
In summary, we calculated the orbital susceptibility for graphene or for an electron system on a two-dimensional honeycomb lattice based on a newly developed general formula. Our result contains all the contributions up to the first order with respect to the overlap integrals between the nearestneighbor atomic orbitals. In contrast to the previous studies, we found the additional contributions of χ occ:1 and χ 2 . In particular, the former gives asymmetry with respect to the sign change of the chemical potential. Furthermore, the physical origin of the new contribution, χ 2 , is clarified in terms of the corrections to the Peierls phase argument. Furthermore, because of the nature of the two-band model, χ inter gives a finite contribution to χ, which is different from the case of the single-band model discussed in the preceding paper. 16) The relative phase θ k between the atomic orbitals on A and B sublattices plays important roles. Interestingly, θ k is directly related to the chirality of the two Dirac cones and also to the integral called the interband Berry connection.
where R = R Ai − R B j and we have used abbreviations such as θ x = ∂ θ k /∂ k x . The terms with the R-summation in (A·6) originate from the integrals between the nearest-neighbor sites, which are in the first order of the overlap integrals. As in the 1s orbital case, we can show that 1 R,pπpπ = x R,pπpπ = y R,pπpπ = 0. 16) As a result, we obtain Eq. (2.24). In the same way, we can show that (µ, ν = x, y) Furthermore, using the general relation which is derived from the k derivative of the equation for u ℓ ′ k in (2.2), we obtain 15) (A·12) From this relation, we have and where we have used the relation (2.24).
Furthermore, from the k µ and k ν derivatives of Eq. (2.2), we obtain From this relation, we can show that 15) and where we have used the relations in (2.24) and (A·10).

Appendix B: Density of states for graphene
In the case of graphene or a honeycomb lattice, the Brillouin zone is also a honeycomb with a size of 8 √ 3π 2 /9a 2 with a being the nearest-neighbor distance, and the system area is L 2 = 3 √ 3a 2 N/2 with N being the total number of unit cells. Therefore, the density of states per area is obtained as dk x dk y δ (±t|γ k |− µ).
(B·5) This is consistent with the fact that the model is equivalent to the free electrons with an effective mass m * = 2h 2 /3ta 2 at the band edge. D(µ) has a diverging peak at µ = ±2t, which corresponds to the van Hove singularity. From the analytical form of (2.21), we obtain