Orbital Magnetism of Bloch Electrons III. Application to Graphene

TL;DR

This paper provides an exact calculation of graphene's orbital susceptibility, revealing additional contributions beyond traditional models and emphasizing the role of interband effects and atomic orbital phases.

Contribution

It applies a rigorous theory to a two-band graphene model, identifying new susceptibility contributions and clarifying the limitations of previous Peierls phase approaches.

Findings

Additional susceptibility contributions from interband and Fermi surface effects.

The phase between atomic orbitals affects orbital susceptibility.

Corrections to the Peierls phase are necessary for accurate modeling.

Abstract

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.

Equations

$$\chi =\frac{e^2}{{\hslash }^2{c}^2}{k}_{\mathrm{B}}T\sum _{\mathbf{k},n}\mathrm{Tr}{\gamma }_x\mathcal{G}{\gamma }_y\mathcal{G}{\gamma }_x\mathcal{G}{\gamma }_y\mathcal{G},$$

$$\chi ={\chi }_{\mathrm{LP}}+{\chi }_{\mathrm{inter}}+{\chi }_{\mathrm{FS}}+{\chi }_{\mathrm{occ}}.$$

$$V(\mathbf{r})=\sum _{{\mathbf{R}}_i}{V}_0(\mathbf{r}-{\mathbf{R}}_i),$$

$$H_{\mathbf{k}}{u}_{\ell \mathbf{k}}(\mathbf{r})={\epsilon }_{\ell }(\mathbf{k})u_{\ell \mathbf{k}}(\mathbf{r}),$$

$$H_{\mathbf{k}}=\frac{{\hslash }^2{k}^2}{2m}-\frac{i{\hslash }^2}{m}\mathbf{k}\cdot \nabla -\frac{{\hslash }^2}{2m}{\nabla }^2+V(\mathbf{r}).$$

$$V_0(\mathbf{r})=-\frac{Z_{\mathrm{eff}}{e}^2}{r},$$

$${\varphi }_{\mathrm{p}\pi }(\mathbf{r})=\frac{1}{\sqrt{24}(a_{\mathrm{B}}^{\ast }{)}^{5/2}}\sqrt{\frac{3}{4\pi }}z{e}^{-r/2a_{\mathrm{B}}^{\ast }},$$

$$a_{\mathrm{B}}^{\ast }=\frac{a_{\mathrm{B}}}{Z_{\mathrm{eff}}}=\frac{a_{\mathrm{B}}}{3.25},$$

$${\Phi }_{\mathrm{p}\pi }(\mathbf{r}-{\mathbf{R}}_i)={\varphi }_{\mathrm{p}\pi }(\mathbf{r}-{\mathbf{R}}_i)-\sum _{j=n.n.}\frac{s}{2}{\varphi }_{\mathrm{p}\pi }(\mathbf{r}-{\mathbf{R}}_j),$$

$$s=\int {\varphi }_{\mathrm{p}\pi }^{\ast }(\mathbf{r}-{\mathbf{R}}_j){\varphi }_{\mathrm{p}\pi }(\mathbf{r}-{\mathbf{R}}_i)d\mathbf{r}.$$

$${\phi }_{\mathrm{A}\mathbf{k}}^{\mathrm{ortho}}(\mathbf{r})=\frac{1}{\sqrt{N}}\sum _{{\mathbf{R}}_{\mathrm{A}i}}{e}^{-i\mathbf{k}(\mathbf{r}-{\mathbf{R}}_{\mathrm{A}i})}{\Phi }_{\mathrm{p}\pi }(\mathbf{r}-{\mathbf{R}}_{\mathrm{A}i}),$$

$${\phi }_{\mathrm{B}\mathbf{k}}^{\mathrm{ortho}}(\mathbf{r})=\frac{1}{\sqrt{N}}\sum _{{\mathbf{R}}_{\mathrm{B}i}}{e}^{-i\mathbf{k}(\mathbf{r}-{\mathbf{R}}_{\mathrm{B}i})}{\Phi }_{\mathrm{p}\pi }(\mathbf{r}-{\mathbf{R}}_{\mathrm{B}i}).$$

$$h_{nm}(\mathbf{k})=\int {\phi }_{n\mathbf{k}}^{\mathrm{ortho}\ast }(\mathbf{r})H_{\mathbf{k}}{\phi }_{m\mathbf{k}}^{\mathrm{ortho}}(\mathbf{r})d\mathbf{r},$$

$$\begin{array}{rl}{h}_{\mathrm{AA}}& =h_{\mathrm{BB}}=E_{\mathrm{p}\pi }+C_{\mathrm{p}\pi \mathrm{p}\pi },\\ h_{\mathrm{AB}}& =h_{\mathrm{BA}}^{\ast }=-t{\gamma }_{\mathbf{k}},\end{array}$$

$$\begin{array}{rl}{C}_{\mathrm{p}\pi \mathrm{p}\pi }& =\int {\varphi }_{\mathrm{p}\pi }^{\ast }(\mathbf{r})\sum _{\mathbf{R}\ne 0}{V}_0(\mathbf{r}-\mathbf{R}){\varphi }_{\mathrm{p}\pi }(\mathbf{r})d\mathbf{r},\\ t& =t_0+s{c}_{\mathrm{p}\pi },\\ t_0& =-\int {\varphi }_{\mathrm{p}\pi }^{\ast }(\mathbf{r}-\mathbf{R})V_0(\mathbf{r}-\mathbf{R}){\varphi }_{\mathrm{p}\pi }(\mathbf{r})d\mathbf{r},\\ c_{\mathrm{p}\pi }& =\int {\varphi }_{\mathrm{p}\pi }^{\ast }(\mathbf{r})V_0(\mathbf{r}-\mathbf{R}){\varphi }_{\mathrm{p}\pi }(\mathbf{r})d\mathbf{r}.\end{array}$$

$${\gamma }_{\mathbf{k}}=\sum _{\mathbf{R}}{e}^{-i\mathbf{k}\cdot \mathbf{R}},$$

$${\gamma }_{\mathbf{k}}=e^{-i{k}_{y}a}+e^{i(\frac{\sqrt{3}}{2}{k}_x+\frac{1}{2}{k}_y)a}+e^{i(-\frac{\sqrt{3}}{2}{k}_x+\frac{1}{2}{k}_y)a},$$

$$\begin{array}{rl}s& =\left(1+\tilde{p}+\frac{2}{5}{\tilde{p}}^2+\frac{{\tilde{p}}^3}{15}\right)e^{-\tilde{p}},\\ t_0& =\frac{Z_{\mathrm{eff}}{e}^2}{4a_{\mathrm{B}}^{\ast }}\left(1+\tilde{p}+\frac{{\tilde{p}}^2}{3}\right)e^{-\tilde{p}},\\ c_{\mathrm{p}\pi }& =-\frac{Z_{\mathrm{eff}}{e}^2}{4a_{\mathrm{B}}^{\ast }}\left\{\frac{2}{\tilde{p}}-\frac{3}{{\tilde{p}}^3}+\left(1+\frac{4}{\tilde{p}}+\frac{6}{{\tilde{p}}^2}+\frac{3}{{\tilde{p}}^3}\right)e^{-2\tilde{p}}\right\},\end{array}$$

$$\begin{array}{rl}{u}_{\mathrm{p}\pi \mathbf{k}}^{+}(\mathbf{r})& =\frac{i}{\sqrt{2}}\left\{e^{\frac{i}{2}{\theta }_{\mathbf{k}}}{\phi }_{\mathrm{A}\mathbf{k}}^{\mathrm{ortho}}(\mathbf{r})-e^{-\frac{i}{2}{\theta }_{\mathbf{k}}}{\phi }_{\mathrm{B}\mathbf{k}}^{\mathrm{ortho}}(\mathbf{r})\right\},\\ & \mathrm{for}{\epsilon }_{\mathrm{p}\pi }^{+}(\mathbf{k})=E_{\mathrm{p}\pi }+C_{\mathrm{p}\pi \mathrm{p}\pi }+{\epsilon }_{\mathbf{k}},\\ u_{\mathrm{p}\pi \mathbf{k}}^{-}(\mathbf{r})& =\frac{1}{\sqrt{2}}\left\{e^{\frac{i}{2}{\theta }_{\mathbf{k}}}{\phi }_{\mathrm{A}\mathbf{k}}^{\mathrm{ortho}}(\mathbf{r})+e^{-\frac{i}{2}{\theta }_{\mathbf{k}}}{\phi }_{\mathrm{B}\mathbf{k}}^{\mathrm{ortho}}(\mathbf{r})\right\},\\ & \mathrm{for}{\epsilon }_{\mathrm{p}\pi }^{-}(\mathbf{k})=E_{\mathrm{p}\pi }+C_{\mathrm{p}\pi \mathrm{p}\pi }-{\epsilon }_{\mathbf{k}},\end{array}$$

$$\begin{array}{rl}{\epsilon }_{\mathbf{k}}& =t|{\gamma }_{\mathbf{k}}|,\\ |{\gamma }_{\mathbf{k}}|& =\sqrt{1+4{\cos }^2\frac{\sqrt{3}{k}_{x}a}{2}+4\cos \frac{\sqrt{3}{k}_{x}a}{2}\cos \frac{3k_{y}a}{2}},\end{array}$$

$$e^{i{\theta }_{\mathbf{k}}}=\frac{{\gamma }_{\mathbf{k}}}{|{\gamma }_{\mathbf{k}}|},(-\pi <{\theta }_{\mathbf{k}}\le \pi ).$$

$$u_{\mathrm{p}\pi \mathbf{k}}^{+}(-\mathbf{r})=u_{\mathrm{p}\pi \mathbf{k}}^{+\ast }(\mathbf{r}),u_{\mathrm{p}\pi \mathbf{k}}^{-}(-\mathbf{r})=u_{\mathrm{p}\pi \mathbf{k}}^{-\ast }(\mathbf{r}).$$

$$\begin{array}{rl}D(\mu )=& \frac{2\sqrt{|\mu |/t}}{3\sqrt{3}{\pi }^{2}t{a}^2}K(\kappa ),\mathrm{for}t<|\mu |<3t,\\ D(\mu )=& \frac{2\sqrt{|\mu |/t}}{3\sqrt{3}{\pi }^{2}t{a}^2\kappa }K(\frac{1}{\kappa }),\mathrm{for}|\mu |<t,\end{array}$$

$$\kappa =\frac{1}{4}\sqrt{\frac{(1+|\mu |/t{)}^3(3-|\mu |/t)}{|\mu |/t}}.$$

$${\gamma }_{\mathbf{k}}\sim \mp \frac{3}{2}|\mathbf{k}-{\mathbf{k}}_0^{\pm }|a{e}^{\pm i\eta },$$

$$\int u_{\mathrm{p}\pi \mathbf{k}}^{\pm †}(\mathbf{r})\frac{\partial u_{\mathrm{p}\pi \mathbf{k}}^{\mp }(\mathbf{r})}{\partial k_{\mu }}d\mathbf{r}=\pm \frac{1}{2}\frac{\partial {\theta }_{\mathbf{k}}}{\partial k_{\mu }}+O(s^2),$$

$${\chi }_{\mathrm{LP}}=\frac{e^2}{6{\hslash }^2{c}^2}\sum _{\mathbf{k},\pm }{f}^{\prime }(\pm {\epsilon }_{\mathbf{k}})\left({\epsilon }_{xx}{\epsilon }_{yy}-{\epsilon }_{xy}^2\right),$$

$${\epsilon }_x=\frac{\partial {\epsilon }_{\mathbf{k}}}{\partial k_x},{\epsilon }_{xx}=\frac{{\partial }^2{\epsilon }_{\mathbf{k}}}{\partial k_x^2},{\epsilon }_{xy}=\frac{{\partial }^2{\epsilon }_{\mathbf{k}}}{\partial k_x\partial k_y},etc.$$

$$\begin{split}\frac{\partial u_{{\rm p}\pi{\bm{k}}}^{\pm}}{\partial k_{x}}&=\frac{C^{\pm}}{\sqrt{2N}}\biggl{[}\sum_{{\bm{R}}_{{\rm A}i}}\left(x-R_{{\rm A}ix}-\frac{1}{2}\theta_{x}\right)\cr&\qquad\qquad\times e^{\frac{i}{2}\theta_{\bm{k}}}e^{-i{\bm{k}}({\bm{r}}-{\bm{R}}_{{\rm A}i})}\Phi_{{\rm p}\pi}({\bm{r}}-{\bm{R}}_{{\rm A}i})\cr&\quad\mp\sum_{{\bm{R}}_{{\rm B}i}}\left(x-R_{{\rm B}ix}+\frac{1}{2}\theta_{x}\right)\cr&\qquad\qquad\times e^{-\frac{i}{2}\theta_{\bm{k}}}e^{-i{\bm{k}}({\bm{r}}-{\bm{R}}_{{\rm B}i})}\Phi_{{\rm p}\pi}({\bm{r}}-{\bm{R}}_{{\rm B}i})\biggr{]},\end{split}$$

$$\begin{split}\chi_{\rm inter}&=-\frac{e^{2}}{\hbar^{2}c^{2}}\sum_{{\bm{k}}}\sum_{\varepsilon_{\ell}=\varepsilon_{{\rm p}\pi}^{+}({\bm{k}}),\varepsilon_{{\rm p}\pi}^{-}({\bm{k}})}\ \sum_{\ell^{\prime}\neq\ell}\frac{f(\varepsilon_{\ell})}{\varepsilon_{\ell}-\varepsilon_{\ell^{\prime}}}\cr&\times\biggl{|}\int\frac{\partial u_{\ell{\bm{k}}}^{\dagger}}{\partial k_{x}}\left(\frac{\partial H_{\bm{k}}}{\partial k_{y}}+\frac{\partial\varepsilon_{\ell}}{\partial k_{y}}\right)u_{\ell^{\prime}{\bm{k}}}d{\bm{r}}-(x\leftrightarrow y)\biggr{|}^{2},\end{split}$$

Full text

1 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 $\mu$ where the two Dirac cones come in contact with each other. This is confirmed by using the exact one-line formula (Fukuyama formula).[7, 8] However, these calculations assume a linear dispersion $\varepsilon=\pm v|{\bm{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 two-dimensional honeycomb lattice was studied[9, 10, 11] using the Fukuyama formula[7]

[TABLE]

where $\cal G$ represents the thermal Green’s function ${\cal G}({\bm{k}},\varepsilon_{n})$ of the matrix form with respect to the band indices, $\varepsilon_{n}$ is the Matsubara frequency, and $\gamma_{\mu}$ ($\mu=x,y$) is the current operator in the $\mu$-direction divided by $e/\hbar$. 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) 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) in the case of the two-band model, the orbital susceptibility for graphene was studied. Calculations based on the Peierls phase in the tight-binding 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) in terms of Bloch wave functions and obtained a new and equivalent formula for the orbital susceptibility $\chi$ as follows:

[TABLE]

The suffixes of $\chi_{\rm LP}$, $\chi_{\rm inter}$, $\chi_{\rm FS}$, and $\chi_{\rm occ}$ denote Landau–Peierls, interband, Fermi surface, and occupied states, respectively.[15] Note that the formula in (2) is exact, as is Eq. (1). As shown in II[16], $\chi_{\rm LP}$ is in the first order with respect to overlap integrals. Thus, we calculate each term in (2) for graphene up to the same order with $\chi_{\rm 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, $\chi_{\rm inter}$ vanishes.[16] In contrast, in graphene, which is a typical two-band model, we show that $\chi_{\rm 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 $\chi_{\rm LP}$. Section 4 is devoted to discussion and summary.

2 Atomic orbital model for graphene

2.1 ${\rm p}_{\pi}$ 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 ${\rm p}_{\pi}$ (or pz) 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, px, and py orbitals forming $\sigma$-bonds. We focus on the contributions from the ${\rm p}_{\pi}$ orbitals of carbon atoms considering that only the ${\rm p}_{\pi}$ band crosses the Fermi energy. ] As in I and II,[15, 16] we assume that the periodic potential $V({\bm{r}})$ is written as

[TABLE]

where ${\bm{R}}_{i}$ represents the positions of carbon atoms forming a honeycomb lattice. Under the periodic potential $V({\bm{r}})$, wave functions are given by $e^{i{\bm{k}}\cdot{\bm{r}}}u_{\ell{\bm{k}}}({\bm{r}})$, where $u_{\ell{\bm{k}}}({\bm{r}})$ satisfies

[TABLE]

with

[TABLE]

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

[TABLE]

where $Z_{\rm eff}$ represents the effective charge, $Z_{\rm eff}=3.25$.[19] The ${\rm p}_{\pi}$ (or ${\rm p}_{z}$) atomic orbital is then given by

[TABLE]

where $a^{*}_{\rm B}$ is the renormalized Bohr radius

[TABLE]

with $a_{\rm B}=\hbar^{2}/me^{2}$.

As in II,[16] we use the orthogonal wave functions[20]

[TABLE]

where the $j$-summation represents the sum over the nearest-neighbor (n.n.) sites of ${\bm{R}}_{i}$ and $s$ is the overlap integral

[TABLE]

Note that $s$ is independent of the direction ${\bm{R}}={\bm{R}}_{j}-{\bm{R}}_{i}$ since the ${\rm p}_{\pi}$ 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, $\phi_{{\rm p}\pi}^{*}({\bm{r}}-{\bm{R}}_{j})\phi_{{\rm p}\pi}({\bm{r}}-{\bm{R}}_{i})$, with ${\bm{R}}_{j}\neq{\bm{R}}_{i}$.

Using these orthogonal wave functions, we consider two linear combinations of atomic orbitals (LCAOs) for the A and B sublattices defined as

[TABLE]

and

[TABLE]

Here $N$ is the total number of sites on each sublattice and ${\bm{R}}_{{\rm A}i}$ (${\bm{R}}_{{\rm B}i}$) represents the position of the site in the A (B) sublattice in the $i$-th unit cell.

We assume that $u_{\ell{\bm{k}}}({\bm{r}})$ can be expanded in terms of $\varphi_{{\rm A}{\bm{k}}}^{\rm ortho}({\bm{r}})$ and $\varphi_{{\rm B}{\bm{k}}}^{\rm ortho}({\bm{r}})$. The mixing of other orbitals is neglected. In this approximation, the matrix elements of $H_{\bm{k}}$ are given by

[TABLE]

where $n,m={\rm A}$ or B. Considering that the nearest-neighbor sites are A and B sublattices, we obtain[16]

[TABLE]

where $E_{{\rm p}\pi}$ is the atomic energy eigenvalue for the pπ orbital and

[TABLE]

The derivations of $t_{0}$ and $c_{{\rm p}\pi}$ in the hopping integral $t$ were discussed and justified in II.[16] $\gamma_{\bm{k}}$ is given by

[TABLE]

where ${\bm{R}}$ are the vectors from a B site to its nearest-neighbor A sites, ${\bm{R}}={\bm{R}}_{{\rm A}j}-{\bm{R}}_{{\rm B}i}$. (Here, we have assumed only the nearest-neighbor hopping integrals, but the extension to longer-range hopping integrals is straightforward.) Explicitly, $\gamma_{\bm{k}}$ becomes

[TABLE]

where $a$ is the distance between the nearest-neighbor sites, i.e., $a=|{\bm{R}}|$.

For the ${\rm p}_{\pi}$-orbital, the integrals can be calculated explicitly as[21]

[TABLE]

with ${\tilde{p}}=a/2a_{\rm B}^{*}$. Note that these integrals do not depend on the direction of ${\bm{R}}$. $s$ and $t$ are in the first order with respect to the overlap integrals, which means that they are proportional to $e^{-{\tilde{p}}}$. Using $a=1.42$A in graphene, we obtain ${\tilde{p}}=4.36$, $s=0.237$, and $t=3.55$eV. The overlap integrals and various matrix elements are shown in Appendix A.

2.2 Energy dispersion and chirality around Dirac points

By diagonalizing the Hamiltonian in (14), we obtain normalized eigenfunctions and energy eigenvalues as follows:

[TABLE]

with

[TABLE]

and

[TABLE]

As we can see from Eq. (19), $u_{{\rm p}\pi{\bm{k}}}^{+}({\bm{r}})$ ($u_{{\rm p}\pi{\bm{k}}}^{-}({\bm{r}})$) is the antibonding (bonding) state between the two sublattices with phase factor $e^{\pm i\theta_{\bm{k}}/2}$. $\varepsilon_{{\rm p}\pi}^{+}({\bm{k}})$ ($\varepsilon_{{\rm p}\pi}^{-}({\bm{k}})$) gives the energy dispersion of the upper (lower) Dirac cone. As in II,[16] the constant energy $E_{{\rm p}\pi}+C_{{\rm p}\pi{\rm p}\pi}$ is included in the chemical potential in the following, and we write the Fermi distribution function $f(\pm\varepsilon_{\bm{k}})$ instead of $f(E_{{\rm p}\pi}+C_{{\rm p}\pi{\rm p}\pi}\pm\varepsilon_{\bm{k}})$ for simplicity.

The phase factors $e^{\pm i\theta_{\bm{k}}/2}$ in $u_{{\rm p}\pi{\bm{k}}}^{\pm}({\bm{r}})$ are determined in order to satisfy

[TABLE]

These relations are required in the case of a centrosymmetric potential and have been used in various steps to derive Eq. (2) in I.[15] To prove Eq. (22), it is necessary to note that each A site at ${\bm{R}}_{{\rm A}i}$ has its partner of the B site satisfying $-{\bm{R}}_{{\rm A}i}={\bm{R}}_{{\rm B}j}$. Using this relation, we can see that $\varphi_{{\rm A}{\bm{k}}}^{\rm ortho}(-{\bm{r}})=\varphi_{{\rm B}{\bm{k}}}^{{\rm ortho}*}({\bm{r}})$ holds, which leads to Eq. (22).

The density of states for the honeycomb lattice was obtained in the context of the phonon density of states.[22] It is given by

[TABLE]

where $K(\kappa)$ is the elliptic integral of the first kind with

[TABLE]

Note that $E_{{\rm p}\pi}+C_{{\rm p}\pi{\rm p}\pi}$ is included in $\mu$. For completeness, its derivation and some limiting cases are given in Appendix B.

Two Dirac points exist at ${\bm{k}}={\bm{k}}_{0}^{\pm}=(\pm\frac{4\pi}{3\sqrt{3}a},0)$, and around these Dirac points, $\gamma_{\bm{k}}$ can be expanded as

[TABLE]

where $\eta$ is the angle between the $k_{x}$ axis and vector $\bm{k}$. Thus, the phase $\theta_{\bm{k}}$ of $\gamma_{\bm{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 $\bm{k}$-derivatives of $\theta_{\bm{k}}$ play important roles since $\theta_{\bm{k}}$ is attached to the wave function $u_{{\rm p}\pi{\bm{k}}}^{\pm}({\bm{r}})$. Furthermore, we find that $\partial\theta_{\bm{k}}/\partial k_{\mu}$ ($\mu=x,y$) is related to an integral as

[TABLE]

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.[25, 23, 24]

3 Orbital susceptibility for Graphene

Using the energy dispersion in the previous section, the Landau–Peierls susceptibility[26, 27] from the ${\rm p}_{\pi}$ orbital is given by

[TABLE]

where we have used the abbreviations

[TABLE]

Note that $\partial\varepsilon^{\pm}_{{\rm p}\pi}({\bm{k}})/\partial k_{\mu}=\pm\partial\varepsilon_{\bm{k}}/\partial k_{\mu}=\pm t\partial|\gamma_{\bm{k}}|/\partial k_{\mu}$, which is in the first order with respect to the overlap integrals. As a result, $\chi_{\rm LP}$ is also in the first order of the overlap integrals. In the following, we calculate $\chi_{\rm inter},\chi_{\rm FS}$, and $\chi_{\rm occ}$ up to the same order.

To evaluate $\chi_{\rm inter}$, $\chi_{\rm FS}$, and $\chi_{\rm occ}$, we use

[TABLE]

with $C^{+}=1$, $C^{-}=-i$ and we have also used abbreviations such as $\theta_{x}=\partial\theta_{\bm{k}}/\partial k_{x}$. In the present model, $\chi_{\rm inter}$ is given by[15, 16]

[TABLE]

where the range of the real-space integral $\int\cdots d{\bm{r}}$ has been extended to the whole system size by using the periodicity of $u_{\ell{\bm{k}}}({\bm{r}})$,[15] and $(x\leftrightarrow y)$ represents terms in which $x$ and $y$ are exchanged. Here, we consider the cases where $\varepsilon_{\ell}=\varepsilon_{{\rm p}\pi}^{+}({\bm{k}})$ and $\varepsilon_{\ell}=\varepsilon_{{\rm p}\pi}^{-}({\bm{k}})$. Because of the Fermi distribution function $f(\varepsilon_{\ell})$, 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, px, and py orbitals.] In contrast, we need to take the summation over $\ell^{\prime}$ ($\ell^{\prime}\neq\ell$) in Eq. (30), 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 (30) vanishes both in the zeroth and first order with respect to the overlap integrals since $\hat{L}_{z}\phi_{\rm 1s}({\bm{r}})=0$. However, in the present two-band model, there is a contribution in the zeroth order that involves the derivatives of the phase $\theta_{\bm{k}}$ even if $\hat{L}_{z}\phi_{{\rm p}_{\pi}}({\bm{r}})=0$. This is in sharp contrast to the single-band case. Since there is a zeroth-order term, the infinite summation of $\ell^{\prime}$ in $\chi_{\rm 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 (30) becomes $\varepsilon_{{\rm p}\pi}^{\pm}({\bm{k}})-\varepsilon_{{\rm p}\pi}^{\mp}({\bm{k}})=\pm 2\varepsilon_{\bm{k}}$ when $\varepsilon_{\ell}=\varepsilon_{{\rm p}\pi}^{\pm}({\bm{k}})$ and $\varepsilon_{\ell^{\prime}}=\varepsilon_{{\rm p}\pi}^{\mp}({\bm{k}})$. 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 $\frac{\partial H_{\bm{k}}}{\partial k_{y}}+\frac{\partial\varepsilon_{\ell}}{\partial k_{y}}$ in Eq. (30) is important.

As shown in Appendix C, the summation over $\ell^{\prime}$ in (30) is carried out analytically, and we obtain

[TABLE]

Note that all the terms are related to the derivatives of $\theta_{\bm{k}}$. In the single-band case,[16] the relative phase $\theta_{\bm{k}}$ does not appear and thus $\chi_{\rm inter}$ vanishes.

Next we calculate $\chi_{\rm FS}$ from the pπ orbital, which is given by[15, 16]

[TABLE]

These integrals can be rewritten with the help of the relation $\hat{L}_{z}\phi_{{\rm p}_{\pi}}({\bm{r}})=0$ in a similar way to $\chi_{\rm inter}$. Then we obtain

[TABLE]

Details of the derivation are shown in Appendix C.

Finally, $\chi_{\rm occ}$ is given by

[TABLE]

[Here again, we do not consider the contributions from the other occupied bands.] In a similar way carried out in our preceding paper,[16] $\chi_{\rm occ}$ becomes

[TABLE]

where Re denotes the real part and the expectation value is defined as

[TABLE]

Note that there is a difference between $\Phi_{{\rm p}\pi}({\bm{r}})$ and $\phi_{{\rm p}\pi}({\bm{r}})$ as in Eq. (9). Using the expectation values in (60), we obtain

[TABLE]

where we have used the definition of $\gamma_{\bm{k}}$ in (16).

We numerically calculate $\chi_{\rm LP}$, $\chi_{\rm inter}$, $\chi_{\rm FS}$, and $\chi_{\rm occ}$ as a function of chemical potential $\mu$ at $T=0$. The results are shown in Fig. 2, in which each contribution is normalized by the Pauli susceptibility $\chi_{0}$ at the band edge ($\mu=\pm 3t$), i.e.,

[TABLE]

where the system size is $L^{2}=3\sqrt{3}a^{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^{*}=2\hbar^{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 $\theta_{x}$ and $\theta_{y}$, which originate from the two-band nature. Furthermore, $\chi_{\rm inter}$ is non-zero, which is in sharp contrast to the single-band case, where $\chi_{\rm 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 $\chi_{\rm inter}$ in Eq. (31) and the first two terms of $\chi_{\rm occ}$ in Eq. (37). However, the first term of $\chi_{\rm inter}$ and the second term of $\chi_{\rm occ}$, both of which diverge as $\mu\rightarrow 0$, exactly cancel with each other. Therefore, only the first term of $\chi_{\rm occ}$ contributes to the total susceptibility in the zeroth order. This term is proportional to the electron number in the pπ band, i.e.,

[TABLE]

where $n(\mu)$ represents the total electron number with spin degeneracy when the chemical potential is $\mu$, and it is calculated from the density of states in (23) as $n(\mu)/L^{2}=2\int_{-3t}^{\mu}D(\mu)d\mu$. $\chi_{\rm occ:1}$ represents the contributions from the occupied states in the partially filled pπ-band, which we call intraband atomic diamagnetism.[16] Since $\chi_{\rm occ:1}$ does not have a factor of $e^{-{\tilde{p}}}$, it gives comparative contributions as $\chi_{\rm LP}$ in a similar way to the single-band case studied in II.

(3) At the band bottom ($\mu=-3t$), only $\chi_{\rm LP}$ has a contribution. Its value is just equal to $-1/3\chi_{0}$, which is understood as Landau’s diamagnetic orbital susceptibility for free electrons. Furthermore, $\chi_{\rm LP}$ has a diverging peak at $\mu=-t$, which corresponds to the van Hove singularity.

(4) $\chi_{\rm FS}$ is always negative and its small wiggle at $\mu=-t$ is owing to a subtle cancellation between the last two terms in (33), both of which diverge as $\mu\rightarrow 0$. $\chi_{\rm FS}$ has a sizable contribution in the region of $-3t<\mu<-t$.

The total susceptibility $\chi=\chi_{\rm LP}+\chi_{\rm inter}+\chi_{\rm FS}+\chi_{\rm occ}$ is shown in Fig. 3 as a function of $\mu$. A $\delta$-function-like peak at $\mu=0$[6] is not included, which we discuss shortly. It is rather surprising that the total of each contribution becomes a smooth function of $\mu$ irrespective of the irregular $\mu$-dependences of $\chi_{\rm inter}$ and $\chi_{\rm FS}$ near $\mu=-t$ in Fig. 2. In Fig. 3, the present result is compared with that obtained by Gòmez-Santos and Stauber,[13] denoted as $\chi^{({\rm GS})}$. Apparently, there is a sizable difference from $\chi^{({\rm GS})}$, which is discussed in the next section.

In the present results, there is asymmetry with respect to the sign change of $\mu$. This is because the contribution $\chi_{{\rm occ}:1}$ in Eq. (39) is a monotonically decreasing function of $\mu$. When the pπ band is fully filled (i.e., $\mu>3t$), only $\chi_{{\rm occ}:1}$ gives the contribution to the orbital susceptibility. Therefore, the total susceptibility becomes

[TABLE]

with $N_{e}$ being the total electron number of the pπ band. This is nothing but the atomic diamagnetism from the pπ electrons.

4 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 $\chi_{\rm FS}$ in (33), we find that the total susceptibility can be rewritten as the following simple form:

[TABLE]

with

[TABLE]

[TABLE]

Note that only $\chi_{\rm occ:1}$ is in the zeroth order with respect to the overlap integrals, while the other three contributions are in the first order. $\chi_{\rm occ:1}$ originates from $\chi_{\rm occ}$, but $\chi_{1}$ and $\chi_{2}$ are combinations of $\chi_{\rm inter}$, $\chi_{\rm FS}$, and $\chi_{\rm occ}$. Finally, when we use the interesting relation

[TABLE]

(see Appendix D), $\chi_{2}$ can be rewritten as

[TABLE]

with

[TABLE]

where we have used the relation in (63). Expectation values for an operator $\cal O$ in terms of $\phi_{{\rm p}\pi}({\bm{r}})$ are defined as

[TABLE]

Using the expectation values calculated in Appendix A, each term in Eq. (46) becomes

[TABLE]

when ${\tilde{p}}=a/2a_{\rm B}^{*}=4.36$.

Each contribution, $\chi_{\rm occ:1}$, $\chi_{\rm LP}$, $\chi_{1}$, and $\chi_{2}$, is plotted in Fig. 4 together with the total susceptibility. The flat part of the total susceptibility near $\mu=0$ originates from the cancellation between the $\mu$-dependences of $\chi_{\rm LP}$ and $\chi_{1}$. As discussed before, only $\chi_{\rm occ:1}$ has asymmetry with respect to $\pm\mu$.

Here we comment on the $\delta$-function-like peak at $\mu=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., $\varepsilon=\pm v|{\bm{k}}|$, we expect

[TABLE]

as a contribution from the Dirac cones at $T=0$. Here $\Gamma$ 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 $\delta$-function-like peak obtained by McClure[6] in clean systems at finite temperatures. This peak will originate from the singular behavior of the $\bm{k}$-derivatives of $\theta_{\bm{k}}$ at the Dirac points. In the present model, we have $v=3ta/2$ from Eq. (25) and thus $\chi^{({\rm Dirac})}$ is proportional to $t^{2}$, which means that $\chi^{({\rm Dirac})}$ will appear in the second order with respect to the overlap integrals. Although $\chi^{({\rm 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:

[TABLE]

where $\cal G$ is now a $2\times 2$ matrix, ${\cal G}=(i\varepsilon_{n}-{\cal H}_{\bm{k}})^{-1}$, and ${\hat{\gamma}}_{\mu}=\partial{\cal H}_{\bm{k}}/\partial k_{\mu}$, etc., with ${\cal H}_{\bm{k}}$ being the Hamiltonian in a $2\times 2$ matrix form. (The spin degeneracy has been included.) The second term in (50) is the “correction term”[13] added to the exact formula of (1).

Actually, before Gòmez-Santos and Stauber, Koshino and Ando[12] used another formula,

[TABLE]

although they did not calculate $\chi^{({\rm KA})}$ explicitly. The last three terms are their “correction terms”. We can show that $\chi^{({\rm GS})}$ and $\chi^{({\rm KA})}$ are equivalent using the relation $\partial{\cal G}_{\bm{k}}/\partial k_{\mu}={\cal G}{\hat{\gamma}}_{\mu}{\cal G}$ and integration by parts.

In the present notation, the $2\times 2$ Hamiltonian is given by

[TABLE]

with $\sigma_{x,y}$ being the Pauli matrix. After some algebra, we can show that $\chi^{({\rm GS})}$ is given in the present notation as

[TABLE]

When we perform integration by parts in the second summation, we can see that $\chi^{({\rm GS})}$ is exactly equal to $\chi_{\rm LP}+\chi_{1}$ in Eq. (41). Therefore, the difference between the present result and $\chi^{({\rm GS})}$ is simply $\chi_{{\rm occ}:1}$ and $\chi_{2}$.

The origin of this difference will be the following. In $\chi^{({\rm GS})}$, (i) the effect of deformation of the wave function $u_{\ell{\bm{k}}}$ is not taken into account completely and (ii) the occupied electron contributions in $\chi_{\rm 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 $\chi^{({\rm GS})}$ and $\chi^{({\rm 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 (4). As a result, the exact Fukuyama formula in Eq. (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 $\chi^{({\rm 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 $\chi_{{\rm occ}:1}$ and $\chi_{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 tight-binding 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 $\Delta t=t_{2}\tilde{h}^{2}$ with $\tilde{h}=eHa^{2}/2c\hbar$.[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

[TABLE]

This formula is equivalent to $\chi_{2}$ in Eq. (45), and $b$ in (45) corresponds to $a^{4}t_{2}/t$. In the present notation, $t_{2}$ is given by[18]

[TABLE]

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}({\bm{r}}-{\bm{R}})$, and we think that this difference is not important. Therefore, $\Delta\chi$ qualitatively corresponds to $\chi_{2}$. To confirm this correspondence, we estimate the second term in (55) as follows:

[TABLE]

where we have used the definition of $t$ in (15). Using the expectation value in (60) and $R_{x}^{2}+R_{y}^{2}=a^{2}$, (56) becomes

[TABLE]

In this estimation, $t_{2}$ becomes

[TABLE]

Substituting this expression into (54), we recover the results in (45) and (46).

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 (2) except for the numerical factor of $\chi_{\rm 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 $\chi_{\rm FS}/2$) and that of the Van Vleck susceptibility (which corresponds to $\chi_{\rm inter}$) vanish in the case of gapped graphene. Furthermore, Langevin-type and geometrical susceptibilities (which correspond to $\chi_{\rm occ}$) give a symmetric contribution with respect to the sign change of $\mu$. 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({\bm{r}})$ is noncentrosymmetric. In such cases, orbital magnetization[29] and the connection to the so-called “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 nearest-neighbor atomic orbitals. In contrast to the previous studies, we found the additional contributions of $\chi_{\rm occ:1}$ and $\chi_{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, $\chi_{2}$, is clarified in terms of the corrections to the Peierls phase argument. Furthermore, because of the nature of the two-band model, $\chi_{\rm inter}$ gives a finite contribution to $\chi$, which is different from the case of the single-band model discussed in the preceding paper.[16] The relative phase $\theta_{\bm{k}}$ between the atomic orbitals on A and B sublattices plays important roles. Interestingly, $\theta_{\bm{k}}$ is directly related to the chirality of the two Dirac cones and also to the integral called the interband Berry connection.

Acknowledgments

We thank H. Fukuyama, F. Piéchon, H. Matsuura, I. Proskurin, T. Kariyado, Y. Fuseya, T. Mizoguchi and N. Okuma for very fruitful discussions. This work was supported by a Grant-in-Aid for Scientific Research on “Multiferroics in Dirac electron materials” (No. 15H02108).

Appendix A Overlap integrals and matrix elements

Various kinds of expectation values defined in eqs. (10), (15), and (47) can be obtained analytically for $\phi_{{\rm p}\pi}({\bm{r}})$. Without loss of generality, we assume ${\bm{R}}=(a,0,0)$. Then using a change of coordinates, $\xi=r+r_{b},\eta=r-r_{b}$ with $r=|{\bm{r}}|,r_{b}=|{\bm{r}}-{\bm{R}}|$,[21] we obtain

[TABLE]

for the pπ orbital with ${\tilde{p}}=a/2a_{\rm B}^{*}$, where $r_{\perp}(r_{\parallel})$ denotes the coordinate in the direction perpendicular (parallel) to $\bm{R}$. In a similar way, we obtain $s$, $t_{0}$ and $c_{{\rm p}\pi}$ in Eq. (18). When we set $\tilde{p}=0$, we have

[TABLE]

Expectation values in terms of $\Phi_{{\rm p}\pi}({\bm{r}})$ defined in Eq. (36) can be easily obtained from $\langle{\cal O}\rangle_{{\rm p}\pi{\rm p}\pi}^{(0)}$. For example,[16]

[TABLE]

and

[TABLE]

Using these relations, we have

[TABLE]

and $\langle x^{2}+y^{2}\rangle_{R,{\rm p}\pi{\rm p}\pi}^{(0)}$ is obtained from $\langle r_{\parallel}^{2}\rangle_{R,{\rm p}\pi{\rm p}\pi}^{(0)}+\langle r_{\perp}^{2}\rangle_{R,{\rm p}\pi{\rm p}\pi}^{(0)}$.

Next we calculate the matrix elements that appear in $\chi_{\rm inter}$, $\chi_{\rm FS}$, and $\chi_{\rm occ}$. For example, using the $\bm{k}$-derivatives of $u_{{\rm p}\pi{\bm{k}}}^{\pm}({\bm{r}})$ in (29), we obtain the interband matrix element as

[TABLE]

where ${\bm{R}}={\bm{R}}_{{\rm A}i}-{\bm{R}}_{{\rm B}j}$ and we have used abbreviations such as $\theta_{x}=\partial\theta_{\bm{k}}/\partial k_{x}$. The terms with the ${\bm{R}}$-summation in (64) 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 $\langle 1\rangle_{R,{\rm p}\pi{\rm p}\pi}=\langle x\rangle_{R,{\rm p}\pi{\rm p}\pi}=\langle y\rangle_{R,{\rm p}\pi{\rm p}\pi}=0$.[16] As a result, we obtain Eq. (26). In the same way, we can show that ($\mu,\nu=x,y$)

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Furthermore, using the general relation

[TABLE]

which is derived from the ${\bm{k}}$ derivative of the equation for $u_{\ell^{\prime}{\bm{k}}}$ in (4), we obtain[15]

[TABLE]

From this relation, we have

[TABLE]

and

[TABLE]

where we have used the relation (26).

Furthermore, from the $k_{\mu}$ and $k_{\nu}$ derivatives of Eq. (4), we obtain

[TABLE]

From this relation, we can show that[15]

[TABLE]

[TABLE]

and

[TABLE]

where we have used the relations in (26) and (68).

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\sqrt{3}\pi^{2}/9a^{2}$ with $a$ being the nearest-neighbor distance, and the system area is $L^{2}=3\sqrt{3}a^{2}N/2$ with $N$ being the total number of unit cells. Therefore, the density of states per area is obtained as

[TABLE]

After performing the $k_{y}$-integral and the change of variable $x=\cos^{2}\sqrt{3}k_{x}a/2$, we obtain

[TABLE]

with

[TABLE]

Finally, using the formula[16, 30]

[TABLE]

with $q=\sqrt{{(a-d)(b-c)}/{(a-c)(b-d)}}$ for $a>b>c>d$, we obtain Eq. (23).

At the bottom (or top) of the band ($\mu=\pm 3t$), $\kappa$ becomes zero and using $K(0)=\pi/2$, we obtain

[TABLE]

This is consistent with the fact that the model is equivalent to the free electrons with an effective mass $m^{*}=2\hbar^{2}/3ta^{2}$ at the band edge.

$D(\mu)$ has a diverging peak at $\mu=\pm 2t$, which corresponds to the van Hove singularity. From the analytical form of (23), we obtain

[TABLE]

where we have used the fact that $\kappa$ behaves as

[TABLE]

near $\mu=1$ and $K(\kappa)\sim\ln(4/\sqrt{1-\kappa^{2}})$ as $\kappa\rightarrow 1$. On the other hand, when $\mu\sim 0$, we can show that

[TABLE]

This is consistent with the density of states of the two massless Dirac electrons around the two Dirac points with velocity $3ta/2$ [see Eq. (25)].

Appendix C Calculation of $\chi_{\rm inter}$ and $\chi_{\rm FS}$ for graphene

First, we calculate the integrals in the second line of (30) for $\chi_{\rm inter}$. Using the $\bm{k}$-derivatives of $u_{{\rm p}\pi{\bm{k}}}^{\pm}$ in (29) and the relation $\hat{L}_{z}\phi_{{\rm p}_{\pi}}({\bm{r}})=0$, we can show that

[TABLE]

with ${\bm{R}}={\bm{R}}_{{\rm A}i}-{\bm{R}}_{{\rm B}j}$. Note that the terms with the $\bm{R}$-summation appear since $\Phi_{{\rm p}\pi}({\bm{r}}-{\bm{R}}_{{\rm A}i})$ contains the nearest-neighbor orbitals, $-(s/2)\phi_{{\rm p}\pi}({\bm{r}}-{\bm{R}}_{{\rm B}j})$ in addition to $\phi_{{\rm p}\pi}({\bm{r}}-{\bm{R}}_{{\rm A}i})$ [see Eq. (9)]. The $\bm{R}$-summation in (85) can be carried out as

[TABLE]

Also, using the definition of $u_{{\rm p}\pi{\bm{k}}}^{\pm}$ in (19), we can rewrite (85) as

[TABLE]

Note that the second term on the right-hand side has $u_{{\rm p}\pi{\bm{k}}}^{\mp}$, which means that the bonding and antibonding orbitals are exchanged in this term. Using this relation, we can see that the matrix element of $\chi_{\rm inter}$ in (30) becomes

[TABLE]

Let us first calculate the case with $u_{\ell^{\prime}{\bm{k}}}=u_{{\rm p}\pi{\bm{k}}}^{\mp}$. In this case, with the help of (71) and (72), the first term in (88) becomes $\theta_{x}\varepsilon_{y}/2+O(s^{2})$, which cancels with the second term using the relation (26). This cancellation means that the numerator in $\chi_{\rm inter}$ in (30) when $u_{\ell^{\prime}{\bm{k}}}=u_{{\rm p}\pi{\bm{k}}}^{\mp}$ is proportional to the fourth order of overlap integrals, i.e., $O(s^{4})$.

Next, we consider the case with $\varepsilon_{\ell^{\prime}}\neq\varepsilon_{{\rm p}\pi\bm{k}}^{\mp}$. In this case, with the help of the general relation in (69), (88) becomes

[TABLE]

where $O(s^{2})$ terms are neglected, and we have used $H_{\bm{k}}u_{{\rm p}\pi{\bm{k}}}^{\pm}=\varepsilon_{{\rm p}\pi}^{\pm}({\bm{k}})u_{{\rm p}\pi{\bm{k}}}^{\pm}$, $\varepsilon_{{\rm p}\pi}^{\pm}({\bm{k}})-\varepsilon_{{\rm p}\pi}^{\mp}({\bm{k}})=\pm 2\varepsilon_{\bm{k}}$, and the relations

[TABLE]

The latter relation is obtained by partial derivative of the former (orthogonality condition).[15]

The factor $(\varepsilon_{{\rm p}\pi\bm{k}}^{\pm}-\varepsilon_{\ell^{\prime}})$ in the first term in (89) cancels with the denominator of $\chi_{\rm inter}$ in (30). As a result, the $\ell^{\prime}$-summation in $\chi_{\rm inter}$ can be carried out using the completeness condition for $u_{\ell^{\prime}{\bm{k}}}$. (A similar method was used in I[15] when obtaining the $f$-sum rule.) Taking account of the fact that the second term in (89) is in the first order with respect to the overlap integrals, we obtain

[TABLE]

where (3 terms) represents the terms in which subscripts $(xxyy)$ are changed to $(xyyx),(yxxy)$, and $(yyxx)$ with minus signs for $(xyyx)$ and $(yxxy)$. Here, we have used (26), (65), (69), (67), (72), (75), and (76). Finally, using (66) and the complex conjugate of (74), we obtain (31).

Next we calculate $\chi_{\rm FS}$. In a similar way to obtain (88), we can show that

[TABLE]

With the help of (74)-(76), (92) becomes

[TABLE]

Here, we have used the transformation in (87) to use the relation in (76) for the last term. Finally, collecting all the terms in $\chi_{\rm FS}$, we obtain (33).

Appendix D Proof of Eq. (44)

Generally we have

[TABLE]

However, $\sum_{{\bm{R}}\neq 0}(R_{x}^{2}+R_{y}^{2}){\rm e}^{-i{\bm{k}}{\bm{R}}}=a^{2}\gamma_{\bm{k}}=a^{2}\varepsilon_{\bm{k}}e^{i\theta_{\bm{k}}}/t$ holds. Therefore, taking the summation of $R_{x}^{2}+R_{y}^{2}$ and taking its real part, we obtain

[TABLE]

On the other hand, its imaginary part gives

[TABLE]

Appendix E Comparison between results of Gómez-Santos et al. and Raoux et al.

The orbital susceptibility obtained by Raoux et al. is given by [see Eq. (30) in Ref.[14].]

[TABLE]

with

[TABLE]

where the subscripts means the partial derivatives with respect to $\bm{k}$. Here, ${\bm{f}}$ is defined as ${\bm{f}}=(t{\rm Re}\ \gamma_{\bm{k}},t{\rm Im}\ \gamma_{\bm{k}},0)$ in the present notation, with Re (Im) representing the real (imaginary) part. Therefore, $|{\bm{f}}|=t|\gamma_{\bm{k}}|$.

In the present notation, $\theta_{\bm{k}}$ is represented as

[TABLE]

Using this relation, we can rewrite $U_{1},U_{2},V_{1}$, and $V_{2}$ as follows:

[TABLE]

Finally substituting them into (97) and using integration by parts, we obtain an orbital susceptibility equal to Gòmez-Santos and Stauber’s result in (53), or equivalently $\chi_{\rm LP}+\chi_{1}$.