Counting Pseudo Landau Levels in Spatially Modulated Dirac Systems
CCounting Pseudo Landau Levels in Spatially Modulated Dirac Systems
Toshikaze Kariyado
International Center for Materials Nanoarchitectonics (WPI-MANA),National Institute for Materials Science, Tsukuba 305-0044, Japan ∗ (Dated: November 12, 2018)In a system with Dirac cones, spatial modulation in material parameters induces a pseudo mag-netic field, which acts like an external magnetic field. Here, we derive a concise formula to count the pseudo Landau levels in the simplest setup for having a pseudo magnetic field. The formula is soconcise that it is helpful in seeing the essence of the phenomenon, and in considering the experimen-tal design for the pseudo magnetic field. Furthermore, it is revealed that anisotropic Dirac cones areadvantageous in pseudo Landau level formation in general. The proposed setup is relatively easyto be realized by spatial modulation in the chemical composition, and we perform an estimationof the pseudo magnetic field in an existing material (an antiperovskite material), by following thecomposition dependence with the help of the ab-initio method.
The external magnetic field is not the only source of theLandau levels. For instance, a certain strain on grapheneleads to the so-called pseudo magnetic field and the resul-tant Landau level structures [1–4]. This phenomenon istied to the most intriguing property of graphene, i.e., theemergent relativistic electron, or the existence of Diraccone in the band structure [5]. Having a Dirac cone, thekey toward the finite pseudo magnetic field is the resem-blance between the shift of the Dirac cone in the Brillouinzone and the minimal coupling in the U(1) gauge the-ory. Since the essence is simply the Dirac cone shift, theidea is not limited to graphene, but is applicable to anysystem with emergent linear dispersion, such as three-dimensional Dirac/Weyl semimetals [6–15].The study of the pseudo magnetic field has several im-portant aspects. Obviously, it is conceptually interestingto observe magnetic phenomena like chiral magnetic ef-fect [12] or quantum oscillations [13] without actuallyapplying magnetic field. Furthermore, the field strengthcan possibly exceed the maximum available strength forthe real magnetic field. As an extreme case, even for asystem inert to the real magnetic field (e.g. neutral par-ticle systems, photonic or phononic crystals [16, 17]), thepseudo magnetic field can be influential as far as thereare Dirac cones. Typically, Dirac cones come with pairs,resulting in multiple Dirac nodes in the Brillouin zone,and the direction of the pseudo magnetic field dependson the nodes. (So, a pseudo magnetic field is regardedas an axial magnetic field.) Then, if the real and pseudomagnetic field coexist, they enhance or cancel with eachother depending on the nodes, which induces the valleyimbalance [13, 18–20]. In that sense, the study of thepseudo magnetic field also has potential importance invalleytronics as next functionalization of materials.As we have noted, the essence of the pseudo magneticfield is the similarity between the Dirac/Weyl node shiftand the minimal coupling. The shift of Dirac/Weyl nodescan be induced in many ways [21–23], and the strain hasbeen frequently used in the context of the pseudo mag-netic field in the literature. However, the strain is not the “ B ”bulk 2bulk 1 bu ff e r l a y e r L xyz k x ∆ k ∆ E bulk 2 bulk 1(a) (b) FIG. 1. (a) (Color online) Setup under consideration. Bulk 1and bulk 2 are similar to each other, but with some relativeshift of the Dirac cones in k x direction. “ B ” represents thedirection of the pseudo magnetic field. (b) Bulk contributionto the band structure as a function of k x . only choice: the spatial dependence of the magnetic mo-ment [6, 24], or the spatially modulating chemical com-position should be equally sufficient. In this paper, weconsider a simple setup for the pseudo magnetic field gen-eration, which is expected to be relatively easy to real-ize with the spatial modulation of the chemical compo-sition. We first give a notably concise formula [Eq. (8)]for the number of observable pseudo Landau levels. Theformula requires only two dimensionless parameters N and R , where N characterizes the length scale of thespatial modulation, while R characterizes the size of theDirac/Weyl node shift. The conciseness of the formulamakes the essence of the pseudo magnetic field transpar-ent, and helps to consider experimental designs. Further-more, inspired by the simple formula, it is pointed outthat anisotropic Dirac cones are better than the isotropicones to appreciate the pseudo Landau level structure. Inthe latter half of this paper, we perform semi ab-initioestimation of R in a real material, antiperovskite A SnO(A=Ca,Sr) [25–28], to make a quantitative argument onthe pseudo magnetic field generation.Let us start with the formulation. The setup in thispaper is illustrated in Fig. 1(a). The considered system a r X i v : . [ c ond - m a t . m e s - h a ll ] J u l consists of three regions, bulk 1, bulk 2, and the bufferregion. For having a transparent discussion, we assumethat bulk 1 (bulk 2) extends to y = + ∞ ( y = −∞ ), whichexcludes free surfaces from our consideration. Bulk 1 andbulk 2 are similar to each other, having Dirac cones in theband structure, but with slightly different node positions.Namely, the effective model for each region is H ( ± ) k = (cid:126) v ( k ± k ) · σ , (1)where k , which denotes the center of Dirac cones, takesdifferent values for the two regions [see Fig. 1(a)]. Inthe buffer region, we assume that the Dirac cones aresmoothly shifted from the position in bulk 1 to the onein bulk 2. For simplicity, our focus is limited to the casethat the Dirac cones are shifted in k x direction. As wecan see from Fig. 1(a), the system is periodic in x and z direction, and we have a band structure as a functionof k x and k z . (If a two-dimensional Dirac/Weyl systemis our target, we simply neglect any structure along z -axis, and omit k z .) Then, the bulk contribution to theband structure typically looks like Fig 1(b) reflecting theshifted node positions. The question is, how about thecontribution from the state at the buffer region.Within the buffer region, the position of the Diraccones depends on y . Then, in a naive treatment, theeffective model is assumed to be H ( ± ) = (cid:126) v ( − i ∇ ± k ( y )) · σ , (2)which is obtained by replacing k by − i ∇ and taking ac-count of the y dependence of k . The latter arguments re-veal that this naive treatment works well. By comparingEq. (2) with the standard minimal coupling ( − i ∇ − e A ),we have A ( ± ) = ∓ (cid:126) e k ( y ) . (3)Once A is given, the pseudo magnetic field is given by B = ∇ × A . In our setup, the node shift is only in k x di-rection, and depends only on y , which means that ∂A x ∂y isthe only relevant component. Assuming that k linearlyinterpolates bulk 1 and bulk 2, the pseudo magnetic fieldstrength B = | B ( ± ) | is estimated as B ∼ (cid:126) e ∆ kL , (4)where L is the thickness of the buffer region, and ∆ k isthe size of the node shift [See Fig. 1(b)].To make the formulation concise, we introduce two di-mensionless parameters N and R respectively as L = N a and ∆ k = 2 πR/a , where a is the lattice constant. (It isimplicitly assumed that the lattice constants in x and y are the same, but the extension to anisotropic casesis trivial.) N represents the length scale of the spatialmodulation, while R measures the size of the node shift in the unit of the Brillouin zone size. As a rough estima-tion, using a typical atomic scale a ∼ B ∼ . × × RN [T] . (5)That is, we potentially have 16 thousand Tesla, and theavailable strength is reduced by a factor of R/N .The obtained pseudo magnetic field induces the Lan-dau levels. Plugging Eq. (4) into a textbook formula, theenergy of n th Landau level is E n = ± (cid:114) πv (cid:126) R | n | N a . (6)[For three-dimensional systems, Eq. (6) corresponds tothe energy at k z = 0.] However, we should note thatnot the all Landau levels are observable in the energyspectrum. That is, since we have the bulk regions aswell as the buffer region, the Landau levels in the bufferregion can be masked by the bulk contribution in theenergy spectrum. It turns out that the diamond regionwith height ∆ E and width ∆ k in Fig. 1(b) is availablefor the Landau levels (see also the latter arguments onthe numerical results in Fig. 2). Since ∆ E is estimatedas ∆ E = (cid:126) v ∆ k , the condition that the n th Landau levelfalls into this diamond region becomes (cid:114) πv (cid:126) R | n | N a < (cid:126) v ∆ k , (7)which leads to a concise expression | n | < π N R. (8)Here we make a short summary: (i) to make the pseudomagnetic field strong, N should be small [Eq. (5)], (ii) toobserve the large number of Landau levels, N should belarge [Eq. (8)], and (iii) large R is always beneficial.Let us move on to the numerical validation of the de-rived formula. For this purpose, we introduce a two-dimensional square lattice tight-binding model with mo-bile Dirac nodes. Specifically, the Hamiltonian is H k = [1 + δ + 2(cos k x + cos k y )] σ z + 2 α sin k y σ y (9)where the lattice constant a is set to 1. By expandingthe Hamiltonian with respect to ˜ k x ≡ k x − π/ k y up to the first order in each of the parameters, we endup with H k ∼ −√ (cid:2) (˜ k x − ˜ δ ) σ z + ˜ αk y σ y (cid:3) , (10)where ˜ δ = δ/ √ α = 2 α/ √
3. That is, ˜ δ behaves as A x , and ˜ α is essentially anisotropy of the Fermi velocity, v y /v x . Therefore, if we assign ˜ δ/ π = 0 .
05 for bulk 1,˜ δ/ π = − .
05 for bulk 2, and the linearly interpolated �������������������� ��� ��� ��� ��� ��� � � � � �� � � � ���� � ������������������ � � � � �� ����������������� � � � � �� ��� ��� ��� ��� ��� ��� � � � ���� � � (a) (b)(c) (d)(e) (f) N = 0 π NR = 0 N = 10 π NR ∼ . N = 30 π NR ∼ . N = 50 π NR ∼ . N = 50 v x /v y ∼ N = 50 v x /v y ∼ . N . For (e) and (f), theanisotropy of the Dirac cone is introduced. value for the buffer region, R = 0 . y direction not to have free edges. Namely, the systemconsists of the repetition of the chunk of bulk 1 – buffer– bulk 2 – buffer.Figures 2(a)–2(d) summarizes the results for ˜ α = 1. Ifthere appears some flat sections in the band structure,they can be regarded as the Landau levels. (A flat sec-tion gives a peak in the density of state.) For N = 0[Fig. 2(a)], i.e., if the change between bulk 1 and bulk 2is sharp, only the zeroth Landau level is identified in theband structure. For N = 1 [Fig. 2(b)], which results in π N R ∼ .
8, we still only see the zeroth Landau level.If N is further increased to π N R ∼ . n = 1 Landau level becomes clearly visible, and we seea small signature of the n = 2 Landau level as well. For N = 50 leading to π N R = 3 . n = 3 is possible.All of these observations confirm the formula Eq. (8).It is worth noting that the expected peak structure inthe (local) density of state should be a key to experimen-tal detection of the pseudo Landau levels. Any measure-ments capable of detecting the density of state, such asSTM/STS as a direct measurement or optical conductiv-ity, might be useful.So far, we have been treating the isotropic Dirac cone.Actually, the anisotropy of the Dirac cone gives signif-icant influence on the observable Landau levels. If the Dirac cone becomes anisotropic, v in the left hand sideof Eq. (7) is replaced by v x v y , and v in the right hand sideby v x . Consequently, the formula Eq. (8) is rewritten as | n | < π v x v y N R. (11)This implies that if we have v x > v y , the number ofobservable Landau levels increases compared with theisotropic case with the same N R . Physically, this is be-cause larger v x means larger ∆ E , and smaller v y meanslarger density of states, both of which is advantageousto observe more Landau levels. The formula Eq. (11)is again confirmed using the toy model by modifying ˜ α .Figures 2(e) and 2(f) shows that for v x /v y ∼ α = 2),the number of the Landau levels is doubled comparingwith the isotropic case, while for v x /v y ∼ . α = 0 . n = 0 and n = 1 Landau levels are clearly seen.Therefore, if one attempts to observe large number ofpseudo Landau levels, it is better to focus on the systemwith anisotropic Dirac cones.Hereafter, we work on the quantitative estimation ofthe pseudo magnetic field in existing materials. Havingformulae Eqs. (8) and (11), the estimation of R is essen-tial, and we derive R in semi ab-initio way. Here, semiab-initio means that we apply the first-principles densityfunctional theory [29] with a small assumption in thecrystal structure to calculate the electronic band struc-ture. Since the two-dimensional cases have been stud-ied in graphene extensively, our focus is on the three-dimensional cases – we take the cubic antiperovskitefamily A SnO (A=Ca,Sr), where the three-dimensionalDirac cones are found on k x , k y , and k z axes in the first-principles calculation. (Note that the experimental stud-ies on this materials are now in progress [30–34].) Thegreat advantage of the antiperovskite family is that we al-ready know the way to tune the electronic structure nearthe Fermi energy at the qualitative level [35]. Namely, itis natural to expect that the Dirac cone shift is realized bypreparing Ca − x ) Sr x SnO and adjusting x [35]. Beforewe proceed, we would like to point out a minor disad-vantage of the antiperovskite family. Strictly speaking,there is a tiny mass gap at the Fermi energy, and there-fore, we have to achieve R such that ∆ E is sufficientlylarger than the mass gap.Here, the electronic structure for 0 < x < x = 0 and x = 1. (The computational de-tails are in parallel with Ref. 35.) Using those results, weconstruct effective models for x = 0 and x = 1, and then,the model for arbitrary x is obtained by interpolating theparameters in the effective model. To be quantitative,the construction of the effective model is conducted usingthe maximally localized Wannier function method imple-mented in Wannier90 package [36]. In practice, we con-struct a 12 orbital model where comes from ( spins) × ( p -orbitals on Sn atom + d x − y / d y − z / d z − x on �������� � � � ��� � ��� � � � ���� � ��� � ���� � ��� � ���� � � � � � ������ � � � � � ������ � � � � �� � � � � � � � � ���� � � � � � � ��� � � � � ��� FIG. 3. (Color online) The band structure ofCa − x ) Sr x SnO obtained with the interpolation usingthe MLWF method. crystallographically equivalent Ca/Sr atoms) = × ( + )[26]. Since we intend to focus on the small variationaround x = 0 .
5, we fix the lattice constant as the averageof the lattice constants for x = 0 and x = 1, which areexperimentally known, throughout the calculation. Theband structures along k x axis for x = 4 / x = 5 / k x ∼ . x = 4 / k x ∼ . x = 5 /
9, where the momentum ismeasured in the unit of 2 π/a . At the end, R is estimatedto be approximately 0.006.As a study complementary to the Wannier orbitalbased interpolation, we perform the calculation with asuperstructure, containing Ca and Sr in a ratio x = 4 / x = 5 /
9. In specific, we consider the superstructurein z direction as shown in the right panels of Fig. 4. Sincethe unit cell is not enlarged in x - y plane, the node shift in k x direction can be discussed in exactly the same footingas the previous paragraph (no Brillouin zone folding in k x and k y direction). Again, assuming that x = 4 / x = 5 / x = 0 .
5, we use the av-eraged lattice constant. The band structures for x = 4 / x = 5 / k x axis are shown in Fig. 4. Fromthis result, we extract the node position for x = 4 / k x = 0 . x = 5 / k x = 0 . R = 0 . ����������������� � � � ���� � ��� � ���� � ��� � ���� � � � � � ����� � � � � ��� � � � � �� � � � � � � � � ���� � ������������������ � � � � � ����� � � � � ��� � � � � �� � � � � � (a)(b) O CaSrSnSrCa
FIG. 4. (Color online) The band structure calculated usingthe enlarged unit cell with the composition Ca − x ) Sr x SnO.The used crystal structures are shown in the right panels. than isotropic ones in appreciating the pseudo Landaulevel structure.The author thanks M. Ogata, H. Takagi and A. Vish-wanath for stimulating discussions. The computation inthis work has been done using the facilities of the Super-computer Center, the Institute for Solid State Physics,the University of Tokyo. This work was supported by aGrant-in-Aid for Scientific Research No. 17K14358. ∗ [email protected][1] F. Guinea, M. I. Katsnelson, and A. K. Geim, Nat Phys , 30 (2010).[2] T. Low and F. Guinea, Nano Lett. , 3551 (2010).[3] N. Levy, S. A. Burke, K. L. Meaker, M. Panlasigui,A. Zettl, F. Guinea, A. H. C. Neto, and M. F. Crommie,Science , 544 (2010).[4] N.-C. Yeh, M.-L. Teague, S. Yeom, B. Standley, R.-P.Wu, D. Boyd, and M. Bockrath, Surface Science ,1649 (2011).[5] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S.Novoselov, and A. K. Geim, Rev. Mod. Phys. , 109(2009).[6] C.-X. Liu, P. Ye, and X.-L. Qi, Phys. Rev. B , 235306(2013).[7] Y. Chen, S. Wu, and A. A. Burkov, Phys. Rev. B ,125105 (2013).[8] M. N. Chernodub, A. Cortijo, A. G. Grushin, K. Land-steiner, and M. A. H. Vozmediano, Phys. Rev. B ,081407 (2014).[9] A. Cortijo, Y. Ferreir´os, K. Landsteiner, and M. A. H. Vozmediano, Phys. Rev. Lett. , 177202 (2015).[10] D. I. Pikulin, A. Chen, and M. Franz, Phys. Rev. X ,041021 (2016).[11] A. G. Grushin, J. W. F. Venderbos, A. Vishwanath, andR. Ilan, Phys. Rev. X , 041046 (2016).[12] A. Cortijo, D. Kharzeev, K. Landsteiner, and M. A. H.Vozmediano, Phys. Rev. B , 241405 (2016).[13] T. Liu, D. I. Pikulin, and M. Franz, Phys. Rev. B ,041201 (2017).[14] E. V. Gorbar, V. A. Miransky, I. A. Shovkovy, and P. O.Sukhachov, Phys. Rev. B , 115422 (2017).[15] S. Tchoumakov, M. Civelli, and M. O. Goerbig, Phys.Rev. B , 125306 (2017).[16] C. Brendel, V. Peano, O. J. Painter, and F. Marquardt,Proc. Natl. Acad. Sci. USA , E3390 (2017).[17] H. Abbaszadeh, A. Souslov, J. Paulose, H. Schomerus,and V. Vitelli, arXiv:1610.06406.[18] A. Chaves, L. Covaci, K. Y. Rakhimov, G. A. Farias, andF. M. Peeters, Phys. Rev. B , 205430 (2010).[19] B. Roy, Z.-X. Hu, and K. Yang, Phys. Rev. B , 121408(2013).[20] M. Settnes, J. H. Garc´ıa, and S. Roche,arXiv:1705.09085.[21] G. Montambaux, F. Pi´echon, J.-N. Fuchs, and M. O.Goerbig, Eur. Phys. J. B , 509 (2009).[22] T. Kariyado and Y. Hatsugai, Phys. Rev. B , 245126(2013).[23] T. Kariyado and Y. Hatsugai, Sci. Rep. , 18107 (2015).[24] Y. Araki, A. Yoshida, and K. Nomura, Phys. Rev. B ,115312 (2016).[25] T. Kariyado and M. Ogata, J. Phys. Soc. Jpn. , 083704(2011). [26] T. Kariyado and M. Ogata, J. Phys. Soc. Jpn. , 064701(2012).[27] M. Klintenberg, J. T. Haraldsen, and A. V. Balatsky,Appl. Phys. Res. , 31 (2014).[28] T. H. Hsieh, J. Liu, and L. Fu, Phys. Rev. B , 081112(2014).[29] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car,C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococ-cioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris,G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis,A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari,F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello,L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P.Seitsonen, A. Smogunov, P. Umari, and R. M. Wentz-covitch, J. Phys.: Condens. Matter , 395502 (2009).[30] Y. F. Lee, F. Wu, R. Kumar, F. Hunte, J. Schwartz, andJ. Narayan, Appl. Phys. Lett. , 112101 (2013).[31] J. Nuss, C. M¨uhle, K. Hayama, V. Abdolazimi, andH. Takagi, Acta Cryst. B , 300 (2015).[32] D. Samal, H. Nakamura, and H. Takagi, APL Materials , 076101 (2016).[33] Y. Okamoto, A. Sakamaki, and K. Takenaka, J. Appl.Phys. , 205106 (2016).[34] M. Oudah, A. Ikeda, J. N. Hausmann, S. Yonezawa,T. Fukumoto, S. Kobayashi, M. Sato, and Y. Maeno,Nature Commun. , 13617 EP (2016).[35] T. Kariyado and M. Ogata, arXiv:1705.08934.[36] A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Van-derbilt, and N. Marzari, Comp. Phys. Commun.178