Diamagnetism of Graphene with Long-Range Scatterers

The diamagnetic susceptibility of disordered monolayer graphene containing scatterers with long-range potential is calculated in a self-consistent Born approximation. Explicit numerical results are obtained for a Gaussian potential. The results show that the delta-function susceptibility in the ideal graphene is broadened by disorder and that the broadening is determined by the condition that states can be mixed with those at the Dirac point. For charged impurities, the susceptibility exhibits a double-peak structure as a function of the Fermi energy at zero temperature, due to the rapid decrease in effective scattering strength with the increase of the Fermi level.


§1. Introduction
The graphene is a monolayer graphite sheet recently fabricated and has been attracting attentions theoretically and experimentally.−10) An important feature is the presence of a topological singularity at k = 0.This leads to the presence of a Landau level at the Dirac point, ε = 0, responsible for the singular diamagnetic susceptibility given by a delta function at the Dirac point. 5,11)The purpose of this paper is to study the diamagnetic response of disordered graphene containing scatterers with long-range potential.
Effects of disorder on the susceptibility have been studied for scatterers with potential range smaller than typical electron wavelength in a self-consistent Born approximation, 12) often used for study of transport properties, 13−16) and within the approximation assuming energy-independent broadening. 17,18)Recently, the scheme based on the self-consistent Born approximation was extended to the case of scatterers with long-range potential. 19)In this paper, we use this new scheme to calculate the diamagnetic susceptibility.We consider scatterers with a Gaussian potential to see explicit dependence on the potential range.It should be mentioned that the orbital magnetism was also studied for related materials, such as graphite intercalation compounds, 20−22) carbon nanotube, 9,23−25) graphene ribbons, 26) few-layer graphenes, 3,12,27) bulk graphite, 28) bismuth, 29−31) and organic compounds. 32)he paper is organized as follows: In §2, following a brief review on the electronic states and the method of calculations of Green's function for long-range scatterers in the self-consistent Born approximation, we discuss the method to calculate the orbital susceptibility.In §3 some examples of numerical results are presented for scatterers with a Gaussian potential.A brief discussion is given in §4 together with results for charged-impurity scattering.A short summary is given in §5.§2.Formulation

Effective-mass description
In a graphene sheet the conduction and valence bands consisting of π orbitals cross at K and K' points of the Brillouin zone, where the Fermi level is located. 33,34)lectronic states of the π-bands near a K point are described by the k•p equation: 2−10) H 0 F (r) = εF (r), (2.1) where F is a two-component wave function, ε is energy, σ = (σ x , σ y ) is the Pauli spin matrix, k = ( kx , ky ) = −i∇ is a wave-vector operator, and γ is a band parameter, given by γ = ( √ 3/2)aγ 0 with a the lattice constant (2.46 Å) and γ 0 the nearest-neighbor hopping integral.
The energy, density of states, and electron concentrations are given by ) respectively, where s = ±1, g v = 2 is the valley degeneracy associated with K and K' points, g s = 2 is the spin degeneracy, ε F is the Fermi energy, and sgn(ε) denotes the sign of ε.

Self-consistent Born approximation
Let us consider Green's function defined by with where • • • means average of impurity configurations, Σ(k, ε) is the self-energy matrix, and v i (r) is the potential of scatterers.We consider scatterers with isotropic potential where v i (q)= v i (q) with q = |q|.
where k = |k| and n is a unit vector in the direction of k, i.e., n = k/k.In a self-consistent Born approximation, the selfenergy matrix becomes where n i is the impurity concentration.The self-consistency equations are , (2.12) where θ is the angle between k and k and X ≡ X(k, ε), X ≡ X(k , ε), etc., for simplicity.

Diamagnetic susceptibility
A theoretical scheme to calculate the diamagnetic susceptibility was developed by Fukuyama. 35)It gives the following expression valid in graphene within the k•p scheme: 17,18,36) with (2.16) where L 2 is the system area and (2.17) Before making explicit evaluation of χ(ε), we consider a sum rule valid for χ(ε) in the present system.Let us consider the integral of χ(ε) over the whole band.Using the fact that F (z) is analytic except along the real axis in the complex plane and decays sufficiently rapidly for large |z|, we have (2.18) where C 1 and C 2 are the contours in the complex z plane shown in Fig. 1.
For the integral along C 2 sufficiently away from the real axis, the integral of χ(ε) is dominated by states with high energies away from the Dirac point ε = 0.In the case of scatterers with long-range potential, scattering effects become negligible for such states apart from those of forward scattering, giving rise to small imaginary part of self-energy independent of energy.Because such forward scattering is unimportant, the integral becomes the same as that of states in the absence of disorder.In ideal graphene, the susceptibility is given by a singular delta function, 11) i.e., Therefore, we have the following sum rule: (2.20) In fact, this sum rule will be confirmed within numerical accuracy by explicit calculations shown in the following.Actually, it is valid even for short-range scatterers. 12)he Feynman diagrams for F (ε) within the selfconsistent Born approximation are shown in Fig. 2. We have (2.21) The subtraction of F (c) (ε) is required for the purpose of avoiding double counting because diagrams represented by F (c) (ε) also appear in F (a) (ε) and F (b) (ε).First, we consider current vertex Ĵx (k, ε) shown in Fig. 2 (c), obtained previously. 19)It satisfies Bethe-Salpeter type equation Ĵx (k, ε) In the right hand side of the above equation, we have terms proportional to σ x , (σ where the coefficients satisfy and U being obtained from U by replacing X and Y with X and Y , respectively.We have further defined Within the similar algebra, the expression for Ĵ0 xy (k, ε) consists of the sum of terms proportional to the followings: where, we have used obtained immediately with the use of σ x σ y = iσ z .We have with , for simplicity.The Bethe-Salpeter type equation (2.27) gives the appearance of terms proportional to those of classes I and II.Further, terms in different classes do not mix with each other.Consequently, we can write the vertex as where K 1 ≡ K 1 (k, ε), etc., for simplicity.The equations for K 1 and K 2 become This shows that K 1 and K 2 are independent of each other.For K 3 , K 4 , etc., we have with and ⎛ ⎜ ⎝ where I 1 ≡ I 1 (k, ε), etc., for simplicity.Then, we have and For diagram F (b) (ε), we have to calculate Ĵyx .The result is given by exchanging σ x and σ y in the expression of Ĵxy and Ĵ0 xy except in (σ•n), etc.As a result, we have immediately (2.41) The vertex functions Ĵxy and Ĵyx are given by a (2,2) matrix.Therefore, they can be represented by four independent parameters rather than six parameters In fact, such a formulation using four parameters is also possible as has been made for bilayer graphene. 37)The over-complete representation used above does not give any undesirable results because of the linearity of the Bethe-Salpeter type equation as long as it is numerically solved by a simple iteration.§3.Numerical Calculations In order to see the explicit dependence on the potential range, we first consider scatterers with a Gaussian potential with potential range d.The scattering strength is char-acterized by dimensionless parameter which is the same as that defined previously. 19)e introduce cutoff energy ε c and wave vector k c through ε c = γk c , where we have roughly ε c ≈ γ 0 , corresponding to the region where the linear dispersion is approximately valid, where γ 0 ∼ 3 eV is the hopping integral between nearest-neighbor carbon atoms within a tight-binding model.It is worth mentioning that the actual value of k c is irrelevant when dk c 1, because only states corresponding to |ε| < ∼ γ/d are mixed with each other and contribute to the susceptibility.This fact will be demonstrated in numerical results presented below.For numerical calculations, we discretize wave vector and numerically solve the self-consistency equation and Bethe-Salpeter type equations iteratively in a manner discussed previously. 19)igure 3 shows calculated (a) diamagnetic susceptibility and (b) density of states as a function of energy for W = 0.15 and for varying potential range d, where unit χ c of the vertical axis is defined by The susceptibility has a sharp peak at the Dirac point, accompanied by a long tail leading to appreciable value up to ε = γ/d denoted by vertical arrows.This feature is most clearly seen in the case of long-range scatterers dk c = 20 and 10, but is actually independent of the potential range as will be demonstrated below.The density of states deviates from that of the ideal graphene without disorder in the same energy region.The integral of χ(ε) over the whole energy agrees with eq.(2.20) within relative errors of a few times 10 −3 , corresponding to inaccuracies of the present numerical calculation.Figure 4 shows results for W = 0.2.The feature of the sharp peak at the Dirac point is less prominent and the long tail behavior has become more prominent.The feature that the susceptibility takes appreciable value and the density of states deviates from that of the ideal graphene without disorder for ε < ∼ γ/d remains the same.
Figure 5 shows results for W = 0.3.The feature of the sharp peak at the Dirac point is completely absent and the susceptibility is now dominated by the long tail behavior.The susceptibility takes appreciable values even for ε > γ/d. Figure 6 shows results for W = 0.5.In this case a small dip appears in the susceptibility shown in Fig. 6(a) at the Dirac point, exhibiting a doublepeak structure and the energy width corresponding to nonvanishing susceptibility is further widened.This widening is closely related to the strong modification of the density of states shown in Fig. 6 (b).
Figure 7 shows the susceptibility in the case of (a) dk c = 10 and (b) 20, respectively, for widely varying W .With the increase of W the sharp peak is weakened and disappears around W = 0.3.Then, the susceptibility has a weak double-peak structure in the narrow range 0.4 ≤ W ≤ 0.5 and becomes essentially flat in the energy range increasing with W for W > 0.5.This feature is exactly the same between dk c = 10 and 20.§4.Discussion In ideal graphene, the diamagnetic susceptibility is given by a delta function, i.e., eq.(2.19).−41) One way is to ascribe it to the presence of zero-energy Landau level at the Dirac point independent of the strength of applied magnetic field.In the presence of scatterers, states at the Dirac point are mixed with those with nonzero energy.In the case of weak scattering, such mixing is possible only for states with wave vector k satisfying k < ∼ 1/d and not for those with k > ∼ 1/d.In dirty cases, states with k > ∼ 1/d start to have some contributions of those at zero energy.This simple picture can well account for the qualitative features of the present numerical results.
The fact that potential-range d is the single relevant parameter mentioned in previous sections is demonstrated by plotting χ(ε) in units of χ d given by as a function of energy measured in units of γ/d in Figs. 8  and 9.This gives a curve independent of dk c .The same universality is valid also for the density of states if being plotted in units of (γ/d)/(2πγ 2 ) = 1/(2πγd).Although not shown here, this scaling starts to become invalid in the case dk c = 2, for which characteristic energy γ/d is already close to cutoff ε c .Figures 8 and 9 contain results (dotted lines) obtained by neglecting all vertex corrections, i.e., by setting Ĵx (k, ε) = σ x , Ĵy (k, ε) = σ y , Ĵxy (k, ε) = σ x Ĝ(k, ε)σ y , and Ĵyx (k, ε)= σ y Ĝ(k, ε)σ x .We may think that this roughly corresponds to the introduction of broadening Γ independent of energy and wave vector.The susceptibility within such an approximation was calculated previously, 17,18) giving ImF (ε+i0 The broadening is given by −ih/2τ , where τ is the relaxation time given by Because of the strong variation of the matrix elements as well as D(ε), the above formula is not directly applicable to the present system.It is true that the main contributions to ImF (ε + i0) come from the region of k corresponding to energy ε = γk and thus the assumption of the constant broadening is nearly valid for the evaluation of ImF (ε+i0).However, the strong variation of 1/τ as a function of ε shows that effective Γ in ImF (ε +i0) varies as a function of ε and the resulting susceptibility significantly deviates from the simple Lorentzian even if vertex corrections are completely neglected.The figures also show that vertex corrections play important roles in the behavior of χ(ε).In fact, the susceptibility without vertex corrections has a sharp peak at the Dirac point with broadening increasing with W and decays rapidly with energy.As a result, the susceptibility becomes negligibly small for ε ∼ γ/d for small W . Vertex corrections modify this feature considerably: In the case of weak disorder W ≤ 0.2, the sharpness of the peak remains roughly the same, although its height is reduced because of the presence of the long tail leading to appreciable magnitude for ε ∼ γ/d.In the case of strong disorder W ≥ 0.3, the sharp peak at the Dirac point disappears and only the long-tail part remains.This significant variation with disorder is obtained only if the vertex corrections are properly taken into account.The sum rule for χ(ε) is satisfied without vertex corrections.
According to the Boltzmann result, the diagonal conductivity is independent of carrier concentration n s as long as the dependence of W on n s is not taken into account.Experimentally, however, the conductivity increases almost linearly with n s for sufficiently large n s , showing that the effective scattering strength in actual graphene on SiO 2 substrate varies considerably with n s .Plausible scatterers giving rise to such strong n s dependence are charged impurities. 42,43)imilar calculations are possible also for charged impurities when the screening effect is included within a Thomas-Fermi approximation.The potential is where κ is the static dielectric constant, which is chosen to be 2.5 in the following, and q s is the Thomas-Fermi screening constant given by at zero temperature.The impurity concentration will be measured in units of which roughly corresponds to the number of unit cells in a unit area for ε c ≈ 3γ 0 , but is by one order of magnitude smaller for present ε c = γ 0 .The scattering strength is a function of the density of states at the Fermi level and therefore should be determined in a self-consistent manner for each Fermi level.
Figure 10 shows some examples of the susceptibility and the density of states as a function of the Fermi energy at zero temperature.The susceptibility exhibits a double-peak structure in the vicinity of the Dirac point.This has nothing to do with the appearance of double peaks for scatterers with Gaussian potential shown above.Figure 11 gives some examples of χ(ε) for different values of the Fermi energies.The variation of the broadening with the Fermi energy causes the double peaks in χ(ε F ) at zero temperature, although each curve for a fixed screened potential does not exhibit double peaks.
The diamagnetic susceptibility of graphene is much larger than the real-spin Pauli paramagnetism determined by the electron mass in vacuum because the g factor is close to two due to very small spin-orbit interaction in carbon systems.−47) §5.Summary In summary, we have calculated the diamagnetic susceptibility of disordered monolayer graphene containing scatterers with long-range potential in the selfconsistent Born approximation.Explicit numerical results have been obtained for the Gaussian potential with range d.The results show that the delta-function susceptibility in the ideal graphene is broadened by disorder and that the broadening is determined by the condition that states can be mixed with those at the Dirac point.In fact, in the weak disorder case the energy width is given by γ/d, but becomes larger in dirty cases.For charged impurities, the susceptibility exhibits a double-peak structure as a function of the Fermi energy at zero temperature due to the rapid decrease in effective scattering strength with the increase of the Fermi level.

Figure Captions
) where θ kk is the angle between k and k and • • • represents the average over k under the condition |k | = |k|.

Fig. 1 (Fig. 2 (Fig. 3 (
Fig. 1 (Color online) Contours C 1 and C 2 in the complex z plane.Contour C 1 is along the real axis just above and below and C 2 is sufficiently away from the real axis.Fig. 2 (a) Diagramatic representation of the diamagnetic susceptibility in the self-consistent Born approximation.(b) The Bethe-Salpeter type equation for vertices Ĵxy and Ĵyx .(c) The Bethe-Salpeter type equation for Ĵx and Ĵy .Fig. 3 (Color online) Calculated (a) diamagnetic susceptibility versus energy and corresponding (b) density of states for W = 0.15.The range of scatterers with a Gaussian potential is chosen as dk c = 2, 5, 10, and 20.The energies corresponding to d/γ are denoted by downward arrows.The dashed line in (b) shows the density of states of ideal graphene.

Fig. 6 (Fig. 7 (
Fig. 6 (Color online) Calculated (a) diamagnetic susceptibility versus energy and corresponding (b) density of states for W = 0.5.Fig. 7 (Color online) Calculated diamagnetic susceptibility versus energy and corresponding density of states for (a) dk c = 10 and (b) 20.The disorder parameter W is varied for wide range.

Fig. 8 (
Fig. 8 (Color online) Calculated (a) diamagnetic susceptibility measured in units of χ d and corresponding (b) density of states in units of 1/(2πγd) are plotted against ε(γ/d) −1 for W = 0.15.Curves for different values of dk c is shifted in the vertical direction by a fixed amount as denoted by horizontal straight lines.The dotted lines in (a) show results obtained by neglecting vertex corrections and the dashed lines in (b) show the density of states of ideal graphene.

Fig. 10 (Fig. 11 (
Fig. 10 (Color online) Calculated (a) susceptibility and (b) density of states as a function of the Fermi energy at zero temperature for scattering by charged impurities with the Thomas-Fermi screening approximation.The screened potential strongly varies with the Fermi energy, leading to a double-peak structure in the vicinity of the Dirac point Fig. 11 (Color online) Calculated energy dependence of the susceptibility for a fixed screened potential at different values of the Fermi energies.n i /n c = 0.005.The susceptibility does not show a doublepeak structure for each value of the Fermi level.

Fig. 1 (
Fig. 1 (Color online) Contours C 1 and C 2 in the complex z plane.Contour C 1 is along the real axis just above and below and C 2 is sufficiently away from the real axis.

Fig. 8 (
Fig. 8 (Color online) Calculated (a) diamagnetic susceptibility measured in units of χ d and corresponding (b) density of states in units of 1/(2πγd) are plotted against ε(γ/d) −1 for W = 0.15.Curves for different values of dk c is shifted in the vertical direction by a fixed amount as denoted by horizontal straight lines.The dotted lines in (a) show results obtained by neglecting vertex corrections and the dashed lines in (b) show the density of states of ideal graphene.