RAPID COMMUNICATIONS PHYSICAL REVIEW B 83, 121201(R) (2011) Characterization of Landau subbands in graphite: A tight-binding study Y. H. Ho, 1 J. Wang, 2 Y. H. Chiu, 2,* M. F. Lin, 2,† and W. P. Su 3,‡ 1 Department of Physics, National Sun Yat-Sen University, Kaohsiung 804, Taiwan 2 Department of Physics, National Cheng Kung University, Tainan 701, Taiwan 3 Department of Physics and Texas Center for Superconductivity, University of Houston, Houston, Texas 77204, USA (Received 2 August 2010; revised manuscript received 1 February 2011; published 8 March 2011) A tight-binding model is diagonalized to elucidate the magnetoelectronic properties of bulk graphite. This method reveals far more significant results than previous treatments based on simplified effective-mass approximation. The wave functions are explicitly presented and their spatial distributions appropriately characterize the Landau subbands. The calculated density of states and optical absorption spectra agree well with experiments. Our accurate results enable one to assess to what degree the bulk graphite resembles monolayer and bilayer graphenes. The inclusion of all interlayer coupling leads to rich spectral features awaiting future experimental verifications. DOI: 10.1103/PhysRevB.83.121201 PACS number(s): 71.15.Ap, 71.70.Di, 78.20.Ls Graphene, a single layer of carbon atoms packed in honeycomb lattice, can be regarded as a component atomic layer of ordinary bulk graphite. The fascinating electronic properties of graphene make it a promising candidate for future nanoelectronics applications, which has stimulated intensive experimental and theoretical studies in the past few years. 1–3 The fabrication of graphene layers and subsequent discovery of their physical properties have also led to renewed interests in bulk graphite. 4–7 In this Rapid Communication, we mainly focus on the Landau level structures of Bernal graphite when subjected to a static magnetic field. In three-dimensional (3D) bulk graphite, a perpendicular magnetic field B 0 ˆ z localizes the planar electron motion into the Landau orbitals, while the motion along the field remains intact. The 3D electronic bands are converted into one-dimensional (1D) ones, the so-called Landau subbands. To evaluate those 1D subbands, previous works are primarily based on the effective-mass approximation that has developed from the Slonczewski Weiss McClure (SWMcC) model. 8–12 This is usually used for the description of Landau states near the Fermi level E F . The accuracy of this simplified model is principally limited to the lowest few Landau subbands. 12–16 Further approximation is made by some authors, Koshino and Ando for example, who have kept only one interlayer coupling γ 1 . 8 With that they were able to decompose the energy spectrum of a multilayer graphene into subcomponents effectively identical to the monolayer or bilayer graphene. Such a decomposition is very appealing and has been fre- quently used in the analysis of experimental data, especially magneto-optical properties. 16,17 However, there is actually little theoretical justification for ignoring other interlayer hopping parameters as evidenced by our accurate results. In addition, the comparison with experiments calls for further improvement in theory. The above considerations have led us to examine solutions of a fairly complete model for graphite, the Peierls tight- binding model, which retains all interlayer couplings. 18–21 The obtained wave functions permit detailed classification and characterization of Landau subbands in terms of their spatial distributions. The variation of the wave function along the band edge K -H is particularly significant. It clarifies the important issue of to what extent graphite behaves effectively like a bilayer or monolayer. Our accurate results also yield other important predictions that await further experimental verification. Overall, we have demonstrated that a good tight- binding model can go a long way toward understanding the Landau subbands. In bulk graphite, the planar hexagonal lattice is described by inequivalent A and B atoms while the interlayer stacking results in four sublattices (A 1 , A 2 , B 1 , B 2 ). At B 0 = 0, the tight-binding Hamiltonian is a 4×4 matrix. Distinct from the SWMcC model, its wave function can be directly expressed as a linear combination of the atomic bases:  = c a |A 1 〉+ c a ′ |A 2 〉+ c b |B 1 〉+ c b ′ |B 2 〉. 22 The associated hopping integrals γ ′ i ’s are converted from the SWMcC parameters γ i ’s, 9 including in-plane hopping integral (〈A 1 |H |B 1 〉= γ ′ 0 =−γ 0 =−3.12 eV), hoppings between the nearest neighboring layers (〈A 1 |H |A 2 〉= γ ′ 1 = γ 1 = 0.38 eV, 〈A 1 |H |B 2 〉= γ ′ 3 = γ 3 = 0.28 eV, 〈B 1 |H |B 2 〉= γ ′ 4 = γ 4 = 0.12 eV), hoppings between the next-nearest neighboring layers (〈B 1 |H |B 1 〉= γ ′ 2 = γ 2 /2 =−10.5 meV, 〈A 1 |H |A 1 〉= γ ′ 5 = γ 5 /2 = −1.5 meV), and the enhancement of site energy on A 1,2 atoms (γ ′ 6 =  − γ 2 + γ 5 = 16 meV). In the presence of a magnetic field, the vector potential induces the Peierls phases to accumulate in the Bloch wave functions. 18–20 Those phases are due to path integrals of vector potential, and introduce an extended periodic boundary condi- tion depending on the magnetic flux φ. A flux quantum φ 0 over a graphite hexagon corresponds to a magnetic field of the order of 79,000 T. For a smaller φ, the unit cell has to contain 8φ 0 /φ (=8R B ) atoms and the Hamiltonian becomes an 8R B ×8R B matrix. For the field reachable in experiments, it requires huge computing time to solve the Hamiltonian. We convert this square matrix to a bandlike one. In this way, computing time can be enormously reduced; even down to 1 T the problem is still tractable. In the case of B 0 = 0, wave function can be expressed as the coefficients on four sublattices: A 1,i , A 2,i , B 1,i , B 2,i (i = 1 ∼ 2R B ). The treatment of bandlike matrices and sublattices can be further applied to other layered systems, other stacking orders, or other external field situations. 7,21,23 Figure 1(a) depicts the k z -dispersed Landau subbands at 60 T between zone boundary points K and H . A series of wave functions distributed among the four constituent sublattices 121201-1 1098-0121/2011/83(12)/121201(4) ©2011 American Physical Society