Nonperturbative Matrix Mechanics Approach to Spin-Split Landau Levels and g-Factor in Spin-Orbit Coupled Solids
aa r X i v : . [ c ond - m a t . m e s - h a ll ] J u l Nonperturbative Matrix Mechanics Approach to Spin-Split Landau Levels andg-Factor in Spin-Orbit Coupled Solids
Yuki Izaki and Yuki Fuseya ∗ Department of Engineering Science, University of Electro-Communications, Chofu, Tokyo 182-8585, Japan (Dated: July 5, 2019)We have proposed a fully quantum approach to non-perturbatively calculate the spin-split Landaulevels and g-factor of various spin-orbit coupled solids, based on the k · p theory in the matrixmechanics representation. The new method considers the detailed band structure and the multibandeffect of spin-orbit coupling irrespective of the magnetic field strength. An application of this methodto PbTe, a typical Dirac electron system, is shown. Contrary to popular belief, it is shown thatthe spin-splitting parameter M , which is the ratio of the Zeeman to cyclotron energy, exhibits aremarkable magnetic-field-dependence. This field-dependence can rectify the existing discrepancybetween experimental and theoretical results. We have also shown that M evaluated from the fandiagram plot is different from that determined as the ratio of the Zeeman to cyclotron energy, whichalso overturns common belief. Change in the Zeeman splitting is the most directobservable consequence of spin-orbit coupling (SOC) insolids. The Zeeman energy, which is the energy differencebetween spin-up and down electrons under a magneticfield B , is usually defined as E Z = gµ B B , where g is theg-factor and µ B is the Bohr magneton. For free electrons, g = 2. On the contrary, for itinerant electrons in solids,it is modified due to the SOC and correlation betweenelectrons [1, 2]. Particularly, large g-factors, such as inZn ( g = 170) [3] and Bi ( g = 1060) [4], are characteristicsof the SOC mechanism with low carrier density. There-fore, measurements of the g-factor or the Zeeman energycan provide rare and valuable information of the SOC insolids.Experimentally, the g-factor can be determined fromquantum oscillations, where the frequency F of an oscil-lation is given by F/B = n + 1 / ± M/ n is theLandau level index.) The spin-splitting parameter M isdefined as the ratio of the Zeeman energy to the cyclotronenergy E c = ¯ hω c . It is usually expressed as M zc = E Z E c = gµ B B ¯ hω c = gm c m e , (1)where m e is the free electron mass, m c is the cyclotronmass and ω c = eB/m c . M zc characterizes the relativeenergy scale of the SOC to the kinetic energy in crystals[5]. When the SOC is negligibly small, it is expected tobe M zc ≪ M zc increases as the impact of the SOC be-comes significant. For M zc = 1, the system is equivalentto the Dirac electrons [6–8]. M zc can also be greater thanunity, depending on the contributions from the higher en-ergy bands [5, 9–11]. In recent years, M has attractedrenewed interest since it is related to Berry’s phase as φ B = πM , which is routinely discussed for topologicalmaterials [12–15].Theoretically, on the other hand, it is extremely chal-lenging to develop a fully quantum framework for thecalculation of the g-factor and M , while considering themultiband effect of SOC. Yafet formulated perturbatively a basic idea of the g-factor [1, 16]. This was followed byseveral investigations on the g-factor, especially for semi-conductors [17–20]. Here, the specific symmetry of eachcrystal was analyzed while perturbatively considering themultiband (more than three bands) effect of SOC. ForLuttinger Hamiltonian, which is the effective Hamilto-nian for the two valence bands of Si or Ge, the spin-splitLandau levels can be computed non-perturbatively usinga specific symmetry operation [1, 21, 22]. Thus far, thenon-perturbative computation of the spin-split Landaulevels has been limited to a two-band model, which isthe minimum model for the interband effect [6, 7]. Themultiband models have been solved using the perturba-tive approach only.Recently, a general analytic formula of g-factor wasobtained on the basis of the relativistic multiband k · p Hamiltonian with perturbation theories (L¨owdin parti-tioning) [5]. This method does not require analysis ofthe crystal, which makes it easily applicable to varioussolids. Although we were able to solve the half-a-centuryold mystery on the large anisotropic g-factor and M zc ofholes in Bi using this method, it is still a perturbativemethod. Further, the g-factor could only be obtained upto an order of O (¯ hω c / ∆) ∼ O ( B ), where ∆ is half of theband gap. This is not sufficient for the cases with strongmagnetic fields and narrow gaps, which have garnered aconsiderable interest recently, particularly for topologi-cal insulators, Weyl, and Dirac fermion systems. In fact,an issue related to the spin-splitting of PbTe, which is atypical narrow gap semiconductor, was reported. Usingthe perturbative g-factor formula, it was predicted that M zc = 0 .
83 [23]. On the contrary, M = 0 . − .
57 experi-mentally, as obtained from the analysis of the Shubnikov-de Haas oscillation [24]. This large discrepancy cannotbe attributed to theoretical error, and points towards afundamental problem with the spin-splitting parameterthat has not been encountered yet.In this Letter, we propose a novel non-perturbativematrix mechanics approach to perform a rigorous calcu- k.p
Hamiltonian π - matrix FIG. 1. Schematic image of the π -matrix method. The kine-matical momentum operator π in the k · p Hamiltonian isexpressed in the matrix mechanics representation. lation of the spin-split Landau levels and the g-factor,regardless of the magnitude of the field and the size ofthe gap. We have termed it the π -matrix method in thisLetter (cf. Fig. 1). It is based on a fully quantum the-ory without any semi-classical assumptions and the Bohr-Sommerfeld quantization rule. Here, we do not need anyspecific analysis of the crystal symmetry, whereas theprevious methods [17–22] need an analysis of the specificsymmetry of each crystal, so that the previous theory isunique for a particular crystal symmetry. This methodcan be easily combined with the band calculations, suchas first principles calculations and tight-binding calcula-tions. Consequently, the detailed band structures and themultiband effect of SOC can be taken into account auto-matically. To test the potential of the π -matrix method,we attempt to resolve a recently raised issue of discrep-ancy in M of PbTe. We show that M exhibits a re-markable dependence on magnetic fields, even though ithas been believed to be invariant against the field [cf.Eq. (1)]. Our result can bridge the existing gap betweentheoretical and experimental studies. Further, we showthat the spin-split parameter derived from the fan dia-gram plot is different from that defined as the ratio ofthe Zeeman to cyclotron energy.It was shown by Luttinger and Kohn [25] that the mo-tion of electrons in a periodic potential and a uniformmagnetic field is described by the following equation, inthe so-called Luttinger-Kohn representation [8, 26, 27]: X ℓ ′ σ ′ (cid:20)(cid:18) ǫ ℓ + π m e (cid:19) δ ℓℓ ′ δ σσ ′ + π · v σσ ′ ℓℓ ′ (cid:21) ψ ℓ ′ σ ′ ( r ) = Eψ ℓσ ( r ) . (2) π = − i ¯ h ∇ + e A is the kinematical momentum operatorunder the magnetic field, where A is the vector potentialand e > r ) can be expanded in terms of the band-edge Blochfunctions u ℓσ ( r ) as Ψ( r ) = P ℓσ ψ ℓσ u ℓσ ( r ), where ℓ and σ indicate the ℓ -th Bloch band and its spin, respectively.It may be noted that σ is not the bare spin, but expressesthe degree of freedom of the Kramers doublet. ǫ ℓ is theband-edge energy of the ℓ -th Bloch band. v σσ ′ ℓℓ ′ is thematrix element of the velocity operator between ψ † ℓσ and ψ ℓ ′ σ ′ . The multiband SOC effect is considered by v σσ ′ ℓℓ ′ in a fully relativistic way, which is the strong merit of k · p theory [25]. In principle, the Hamiltonian in the Luttinger-Kohn representation can be written in a matrixform. For a L band system, it is expressed in terms ofa 2 L × L matrix. Initially, it seems that the energy ofsuch electrons can be easily obtained by diagonalizing theHamiltonian, which is true for the cases in the absence ofmagnetic field. However, this is not true in the presenceof magnetic field since the commutation relation π × π = − ie ¯ h B (3)cannot be satisfied by the simple diagonalization. Thiscommutation relation prevents the rigorous calculationof the Landau levels. As discussed earlier, a possible wayto circumvent this theoretical difficulty is the inclusionof the perturbation theory [1, 5, 16–20, 23].We use an unconventional yet simple idea to resolvethis difficulty, which is the incorporation of matrix me-chanics [28, 29]. Because the commutation relation Eq.(3) is essentially the same as that of the harmonic oscil-lation, π can be expressed in terms of a matrix [29]. Byreplacing π in Eq. (2) by the matrix form ˆ π , we obtainthe Hamiltonian as the combination of the Luttinger-Kohn and the matrix mechanics representation (see Fig.1). The basis of the matrix ˆ π is the Landau level index n in the present case. With this matrix of ˆ π , the spin-splitLandau levels can be calculated rigorously by a simplenumerical diagonalization, even though the matrix sizeof the Hamiltonian becomes large, i.e. (2 LN ) × (2 LN ),where N is the number of Landau levels considered. Thismethod is gauge invariant, which guarantees the validityof the theory under a magnetic field, since the commu-tation relation Eq. (3) holds for arbitrary gauge. Theparameters ǫ ℓ and v σσ ′ ℓℓ ′ in the Hamiltonian can be di-rectly calculated from the band calculations, such astight-binding and first principles calculations, at zerofields . To the best of our knowledge, this is the first non-perturbative method that can be used to calculate thespin-split Landau levels while considering the detailedband structure and the multiband effect of SOC underany magnetic field strength. We verified the validity ofthe present method for the well-known case of free elec-trons and Dirac electrons. The details of the π -matrixmethod are given in the Supplemental Material [30].We now apply the π -matrix method to PbTe, which is anarrow gap IV-VI semiconductor with a strong SOC. Asmentioned earlier, an issue regarding large discrepancybetween experimentally and theoretically reported valuesof M in PbTe has been raised recently. To resolve thisissue, we calculate the Zeeman energy, the cyclotron en-ergy, and their ratio M zc , by the π -matrix method alongwith the relativistic tight-binding model by Lent et al [31]. The results obtained agree with those obtained byanother model by Lach-hab et al . [32]. The details aregiven in [30].Figure 2 (a) shows the spin-split Landau levels for thetop valence band at the L -point as a function of the mag-netic field. (The experimental values of M were obtained -160-140-120-100-80-60-40-200 E ( m e V ) E ( m e V ) B (T) M zc π-matrixex-Dirac CyclotronZeemanπ-matrixex-Dirac π-matrix 0-0+1-1+2-2+3-3+(a)(b)(c) ex-Dirac FIG. 2. Magnetic-field-dependence of (a) spin-split Landaulevels, (b) Zeeman E Z and cyclotron E c energy, and (c) spin-splitting parameter M zc . ( E Z,c and M zc are obtained for n = 2.) The solid and dashed lines represent the resultsobtained by the π -matrix method and the extended Diracmodel, respectively. for hole doped samples [24].) The magnetic field is alongthe (001) direction, and the wavenumber along the fieldis set to be at the extremum of energy. The Landau levelindex n and their spins σ are uniquely determined fromthe weak-field-limit value of M zc obtained by using theL¨owdin partitioning [30]. Each energy level exhibits asublinear field-dependence, suggesting that the system isclose to the Dirac electrons, whose energy is given by E D j = p ∆ + 2∆¯ hω c j , ( j = n + σ/ , , , · · · ). Forperfect Dirac electrons, the lowest Landau level is fieldinvariant. However, the obtained lowest Landau level ex-hibits a weak field-dependence. This deviation from theperfect Dirac electrons is a clear indication of the con-tribution from the other bands, i.e., the multiband effectof SOC, which are ignored in the ordinary Dirac electron model. This multiband effect of SOC can be expressed byintroducing the additional g-factor term up to the orderof O ( B ) in the form E exD n,σ = p ∆ + 2∆¯ hω c ( n + σ/
2) + σ g ′ µ B B, (4)which is called as the extended Dirac model [4, 33–35].(For holes, a negative sign is needed in Eq. (4).) Todate, the value of g ′ has been phenomenologically deter-mined to fit the experimental data [4, 35, 36]. We, forthe first time, determine this value microscopically. Thecorresponding results for g ′ = − . n, σ ) = (0 , − ) increases withthe magnetic field when g ′ < B <
60 T.From the spin-split Landau levels, we evaluate the cy-clotron energy E c,n = E n, + − E n − , + and the Zeemanenergy E Z,n = E n, + − E n, − , which are shown in Fig. 2(b) for n = 2. Both E c,n and E Z,n are proportional to B in the weak field limit ( B < ∼ g = 41 .
0, which agrees with the previous valueobtained perturbatively [23, 30]. However, E Z,c exhibita sublinear behavior in the strong field regime. There-fore, the g-factor, which is a coefficient of B , cannotbe well-defined in this regime. This field-dependence isroughly explained by the extended Dirac model, as shownin Fig. 2 (b). The observed deviation of the extendedDirac model from the π -matrix method, which is partic-ularly noticeable for B > ∼
30 T, indicates that the formeris not adequate for this field region. These deviationsarise because the higher order corrections, which are notconsidered in the extended Dirac model but rigorouslyconsidered in the π -matrix method, are not negligible athigh fields. Next, we calculate the spin-splitting param-eter M zc = E Z /E c to investigate the long established as-sumption that M zc is independent of the magnetic field.This is because, in Eq. (1), both E c and E Z have beenobtained up to only O ( B ) in the existing theories. Con-sequently, the field-dependence of M zc has never beenexamined in detail.However, we observe that M zc exhibits a remarkablefield-dependence, which is clear from Fig. 2 (c), where n = 2. At zero field, M zc = 0 .
80, and it rapidly decreasesas M zc = 0 .
46 at B = 55 T. This drastic reduction isremarkable, and we found that it can be qualitativelyexplained using the extended Dirac model. From thismodel, it is easy to derive that M zc = E D j +1 − E D j + g ′ µ B BE D j +1 − E D j . (5)The above equation clearly shows that M zc < B for g ′ <
0. Therefore, the B (T) M z c x=0.00x=0.20x=0.40x=0.60x=0.80x=1.00 FIG. 3. Spin-splitting parameter M zc for Pb − x Sn x Te. Theband inversion (therefore the topological transition) occurs at x = 0 . Landau level index / B n ( / T ) -200 -100 0 μ (meV) M f an FIG. 4. Plot of the reciprocal fields 1 /B n for PbTe, where µ touches the bottom of the n -th Landau level. Here M fan = x − − x + = 0 .
63. The inset shows the µ -dependence of M fan . reduction of M zc in Fig. 2 (c) is a logical consequence ofthe increase in the lowest Landau level in Fig. 2 (a).It is expected that M zc should be an increasing func-tion of B for g ′ >
0. This can be verified by substitut-ing Sn for Pb, i.e., Pb − x Sn x Te. It is well known thatthe band inversion between the conduction and valencebands occurs at around x ≃ . et al. model, the inversion point is at x = 0 .
38, where M zc changes from M zc < M zc > M zc forPb − x Sn x Te using the Lent et al. model. As expected, M zc is an increasing function of B after the inversionat x = 0 .
38, since M zc > g ′ > x > . M is the fan dia-gram plot [2], where 1 /B n is plotted as a function of theLandau level index. Here, B n is the field at which thechemical potential touches the n -th Landau level. These reciprocal fields are linearly fitted for each spin, followingwhich M is evaluated from the difference of x -intercept,i.e., M fan = x − − x + . In principle, M fan is not exactlyequal to M zc . The significant difference can be easily un-derstood from the fact that M fan is independent of B ,but dependent on µ , and vice versa is true for M zc . Therelation M fan = M zc is only true for the cases of freeelectrons with Zeeman energy and Dirac electrons [30].Figure 4 shows the fan diagram of PbTe for the Lent etal. model with µ = −
92 meV from the top of the va-lence band (the hole density is n h = 3 . × cm − ).We obtain M Lentfan = 0 .
63 by fitting the calculated 1 /B n .The µ -dependence of M fan is shown in the inset of Fig. 4. M fan tends to approach the zero-field value of M zc = 0 . µ = 0 limit. This is also consistent with the anal-ysis using the extended Dirac model [30].Experimentally, M zc can be determined directly fromthe Zeeman and cyclotron energy [5]. It is not appro-priate to compare the approximate value of M obtainedfrom the Lifshitz-Kosevich formula with the present the-oretical value, since this formula is not based on a fullyquantum theory. According to Akiba et al. [24], M expfan =0 .
56 (evaluated from Fig. 9), which is less than thatobtained using the Lent et al. model M Lentfan = 0 .
63. Al-though the agreement is not perfect yet, we have success-fully removed the large discrepancy between experimentand theory. Further enhancements of the theoretical ac-curacy can be achieved by improving the accuracy of theband calculations, since the value of M is quite sensitiveto the details of the band structure. Even the high-energybands (more than 1 eV far from the conduction band) canchange M from M Lentfan = 0 .
63 to M Lach − habfan = 0 .
40 [30].In summary, we have demonstrated a novel non-perturbative method, which is based on matrix mechan-ics, for calculating the spin-split Landau levels. This isthe first method to elucidate the following properties ofthe spin-splitting parameter M . (i) M zc is largely de-pendent on the magnetic field; (ii) In general, M zc is notequivalent to M fan . The origin of these previously un-known properties is the multiband effect of SOC and thehigher order corrections in B . This method also pro-vides an explanation for the large discrepancy betweenexperimentally and theoretically reported values of M in PbTe. Apart from the specific case of PbTe consid-ered here, this method can be useful for other systems inwhich the SOC plays a relevant role, such as topologicalinsulators ( M ∼ Se has not been explained yet[9–11]), Weyl, and Dirac fermion systems.We thank M. Tokunaga and K. Akiba for helpful dis-cussions. This work is supported by JSPS KAKENHI(grants No. 16K05437 and 19H01850). ∗ [email protected] [1] Y. Yafet, Solid State Phys. , 1 (1963).[2] D. Shoenberg, Magnetic Oscillations in Metals , Cam-bridge Monographs on Physics (Cambridge UniversityPress, 1984).[3] W. J. O’Sullivan and J. E. Schirber,Phys. Rev. , 519 (1967).[4] Z. Zhu, B. Fauqu´e, Y. Fuseya, and K. Behnia, Phys.Rev. B , 115137 (2011).[5] Y. Fuseya, Z. Zhu, B. Fauqu´e, W. Kang, B. Lenoir, andK. Behnia, Phys. Rev. Lett. , 216401 (2015).[6] M. H. Cohen and E. I. Blount, Phil. Mag. , 115 (1960).[7] P. A. Wolff, J. Phys. Chem. Solids , 1057 (1964).[8] Y. Fuseya, M. Ogata, and H. Fukuyama, J. Phys. Soc.Jpn. , 012001 (2015).[9] H. Kohler and E. Wuchener, Phys. Stat. Sol. , 665(1975).[10] B. Fauqu´e, N. P. Butch, P. Syers, J. Paglione, S. Wied-mann, A. Collaudin, B. Grena, U. Zeitler, andK. Behnia, Phys. Rev. B , 035133 (2013).[11] M. Orlita, B. A. Piot, G. Martinez, N. K. S. Kumar,C. Faugeras, M. Potemski, C. Michel, E. M. Hankiewicz,T. Brauner, C. Draˇsar, S. Schreyeck, S. Grauer, K. Brun-ner, C. Gould, C. Br¨une, and L. W. Molenkamp,Phys. Rev. Lett. , 186401 (2015).[12] D. Xiao, M.-C. Chang, and Q. Niu,Rev. Mod. Phys. , 1959 (2010).[13] Y. Ando, J. Phys. Soc. Jpn. , 102001 (2013).[14] H. Murakawa, M. S. Bahramy, M. Tokunaga, Y. Kohama,C. Bell, Y. Kaneko, N. Nagaosa, H. Y. Hwang, andY. Tokura, Science , 1490 (2013).[15] T. Liang, S. Kushwaha, J. Kim, Q. Gibson,J. Lin, N. Kioussis, R. J. Cava, and N. P. Ong,Science Advances (2017), 10.1126/sciadv.1602510.[16] Y. Yafet, Phys. Rev. , 679 (1957).[17] L. M. Roth, B. Lax, and S. Zwerdling, Phys. Rev. ,90 (1959).[18] M. Cardona, Journal of Physics and Chemistry of Solids , 1543 (1963).[19] C. R. Pidgeon and R. N. Brown,Phys. Rev. , 575 (1966).[20] D. L. Mitchell and R. F. Wallis, Phys. Rev. , 581 (1966).[21] J. M. Luttinger, Phys. Rev. , 1030 (1956).[22] R. F. Wallis and H. J. Bowlden,Phys. Rev. , 456 (1960).[23] H. Hayasaka and Y. Fuseya,Journal of Physics: Condensed Matter , 31LT01 (2016).[24] K. Akiba, A. Miyake, H. Sakai, K. Katayama,H. Murakawa, N. Hanasaki, S. Takaoka,Y. Nakanishi, M. Yoshizawa, and M. Tokunaga,Phys. Rev. B , 115144 (2018).[25] J. M. Luttinger and W. Kohn, Phys. Rev. , 869 (1955).[26] R. Winkler, Spin-Orbit Coupling Effects in Two-Dimensional Electron and Hole Systems (Springer-Verlag, 2003).[27] L. C. L. Y. Voon and M. Willatzen,
The k · p Method (Springer-Verlag, 2009).[28] W. Heisenberg, Zeitschrift f¨ur Physik , 879 (1925).[29] M. Born and P. Jordan,Zeitschrift f¨ur Physik , 858 (1925).[30] See Supplemental Material [url] for details.[31] C. S. Lent, M. A. Bowen, J. D. Dow, and R. S. Allgaier,Superlattices Microstruct. , 491 (1986).[32] M. Lach-hab, M. Keegan, D. Pa-paconstantopoulos, and M. Mehl,Journal of Physics and Chemistry of Solids , 1639 (2000).[33] G. A. Baraff, Phys. Rev. , A842 (1965).[34] M. P. Vecchi, J. R. Pereira, and M. S. Dresselhaus, Phys.Rev. B , 298 (1976).[35] Z. Zhu, B. Fauqu´e, K. Behnia, and Y. Fuseya,J. Phys.: Condens. Matter , 313001 (2018).[36] Z. Zhu, B. Fauqu´e, L. Malone, A. B.Antunes, Y. Fuseya, and K. Behnia,Proc. Natl. Acad. Sci. U.S.A. , 14813 (2012).[37] J. O. Dimmock, in The Physics of Semimetals andNarrow-Gap Semiconductors , edited by D. L. Carter andR. T. Bate (Pargamon, Oxford, 1971) pp. 310–330.[38] T. H. Hsieh, H. Lin, J. Liu, W. Duan, A. Bansil, andL. Fu, Nat Commun3