Table Of Content

Topological Phases in a Weyl Semimetal Multilayer Kazuki Yokomizo1 and Shuichi Murakami1,2 1Department of Physics, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo, 152-8551, Japan 2TIES, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo, 152-8551, Japan We investigate multilayers of a normal insulator and a Weyl semimetal using two models: an effective model and a lattice model. As a result, we find that the behavior of the multilayers 7 is qualitatively different depending on the stacking direction relative to the displacement vector 1 between the Weyl nodes. When the stacking direction is perpendicular to the displacement vector 0 betweentheWeylnodes,thesystemshoweitherthenormalinsulatorortheWeylsemimetalphases 2 dependingonthethicknessesofthetwolayers. Incontact,whenthestackingdirectionisparallelto n that, the phase diagram is rich, containing the normal insulator phase, the Weyl semimetal phase a andthequantumanomalousHallphaseswithvariousvaluesoftheChernnumber. Asasuperlattice J period increases, the Chern number in the quantum anomalous Hall phases increases. Thus, one 9 candesignaquantumanomalousHallsystemwithvariousChernnumbersinamultilayerofaWeyl 1 semimetal and a normal insulator. Applications to Weylsemimetal materials are discussed. ] l l I. INTRODUCTION displacementvectorbetweenthe Weyl nodes. Therefore, a wefocus ontwopatternsofstackingdirections,Patterns h - Recently,varioustopologicalsystemsareactivelystud- A and B shown in Fig. 1 (b) and (c), and calculate the s ied in the field of condensed matter physics. After the- phase diagrams of Patterns A and B separately. In the e m oretical proposals of the quantum spin Hall system1–3, calculation,weuseaneffectivemodelandalatticemodel, various topological systems have been theoretically pre- and their results well agree. In the multilayer with Pat- . t dicted or experimentally observed in real materials. In tern A, the phase diagram have only two phases, a NI a m particular, a Weyl semimetal (WSM) which belongs to and a WSM, depending on the ratio of the thicknesses a class of topological semimetals has been proposed4,5. of the two layers. On the other hand, in the multilayer - d A WSM has Dirac cones without spin degeneracy at withPatterB,the phasediagramcontainsaNI, a WSM n or near the Fermi energy. It cannot be realized un- and QAH phases, depending on the thicknesses of the o less the time-reversal symmetry or the inversion sym- two constituent layers. The Chern number of the QAH c metry is broken. The vertex of the Dirac-cone disper- phases becomes larger and larger when the superlattice [ sion is called a Weyl node, and it can be regarded as periodicity becomes larger. This result offers us a way 1 a monopole or an anti-monopole for the Berry curva- to design QAH systems in a WSM-NI multilayer in a v ture field in the wavevector space4,6,7. Various WSMs controlled way. 2 are predicted after these works, and as candidates of This paper is organized as follows. In Sec. II, we cal- 2 WSMs, pyrochlore iridates5,8,9, HgCr Se 10,11, TaAs12, culate phase diagrams of multilayers for Patterns A and 3 2 4 5 Co2TiX (X=Si,Ge,Sn)13, ZrCo2Sn14, MoxW1−xTe215, Bfromthe effectivemodelforaWSM.We alsocalculate 0 Ge2Sb2Te516, TlBiSe217, Hg1−x−yCdxMnyTe18 in a them using the lattice model and compare the results . magneticfield,andBi Sb inamagneticfield19have of the effective model with those of the lattice model in 1 0.97 0.03 beenproposed. Inrecentyears,someWSMshasbeenex- Sec. III. Finally, we discuss the results and applications 0 7 perimentally observed, such as TaAs20–24 and NbAs25. to real materials in Sec. IV. 1 Among theoretical works of the WSM, a multilayer : of a topological insulator (TI) and a normal insulator v II. MULTILAYER FROM THE EFFECTIVE (NI) with magnetization has been proposed to realize Xi the WSM phase26. Moreover, a TI-NI multilayer with MODEL r broken inversion symmetry by an electric field has also a A. Effective model for a WSM been proposed to show the WSM phase by tuning some parameters,forexamplethethicknessofthelayersinthe multilayerofHgTe/CdTe27. Inaddition,physicalproper- In our calculation on multilayers, we use the effective tiesofmultilayersofaTIandaNIhavebeeninvestigated model for a WSM proposed in Ref. 31. The model de- theoretically28–30. scribes both a WSM and a NI by changing a single pa- In this paper, we study multilayers of a WSM and a rameter m. The Hamiltonian of the model is a 2 2 NI, which has not been studied thus far, to the authors’ matrix H(k,m), which acts on a space consisting o×f a knowledge. We find that the resulting phase diagrams single conduction band and a single valence band. The showrichphysics,includingthequantumanomalousHall 2 2 Hamiltonian H(k,m) is given by × (QAH) phases with various values of the Chern number. H(k,m)=γ k2 m σ +v(k σ +k σ ), (1) Thephasediagramsarequalitativelydifferentdepending x− x y y z z onthestackingdirectionofthemultilayersrelativetothe whereσi (i=1,2,3(cid:0))areth(cid:1)ePaulimatrices,andv andγ 2 arenonzeroconstants. Wechoosethemtobepositivefor (a), and the Hamiltonian is given by simplicity. Energyeigenvaluesofthesystemaregivenby H =γ k2 m(z) σ +v(k σ i∂ σ ), (3) x− x y y − z z E = γ2(k2 m)2+v2 k2+k2 . (2) m(z)=(cid:0) m+ ((cid:1)0≤z ≤a), (4) ± x− y z ( m− ( b z 0), q − − ≤ ≤ (cid:0) (cid:1) m(z+(a+b))=m(z), (5) WeassumethattheFermienergyisatE =0. Whenm< 0, the system is an insulator with the bulk gap 2γ m. wherem± arepositiveconstants. We assumethe system On the other hand, when m > 0, the system is a WS|M|. to be infinite along the xy plane. Along the z axis the The bulk gapcloses attwoWeyl nodes k=( √m,0,0), systemhasperiodicitya+b. Thismodelissimilartothe whichareamonopoleandanantimonopolefo±rtheBerry Krönig-Penney model, and it can be solved similarly, by curvature field. constructing plane-wave solutions in each layer and im- posing proper boundary conditions at the interfaces. In SurfaceFermiarcsofaWSMarealsodescribedwithin the solution, a Bloch wavenumber k along the z direc- this model31. In order to show the surface states, we z tionis introduced. Detailedcalculationsaresummarized consider two types of surfaces of a slab; a top surface in Appendix A1. and a bottom surface, perpendicular to the z direction. We first investigate phases in this multilayer. When The top surface represents the z = 0 surface where the a/b>m /m , the gap at E =E =0 closes at z < 0 region is a WSM with m > 0 and the z > 0 − + F region is the vacuum with m<0. Similarly, the bottom surface represents the z = 0 surface where the z > 0 m a m b region is a WSM and the z < 0 region is the vacuum. k= + − − ,0,0 , (6) As shown in Ref. 31, surface states exist both on the ±r a+b ! top and the bottom surfaces within √m < k < √m. x − These surface states form Fermi arcs. In particular, at andtheyareshowntobetheWeylnodes. Therefore,the E = 0, the surface states at each surface are located at multilayer is in the WSM phase. When a/b < m−/m+, √m < k < √m, k = k = 0, and they connect the the system is gapped. It is the NI phase because it is x y z − Weyl nodes projected on the xy plane. connected to the NI phase when a 0. By increasing → a across the phase boundary a/b = m /m , the Weyl As mentioned earlier, a WSM is realized only when − + nodes are pairwise created from the origin k = 0 in the either the time-reversal symmetry or the inversion sym- wavevector space. By varying the parameters a, b, m , metry is broken. Eq. (1) itself preserves the inversion + or m , the Weyl nodes move along the k axis from the symmetrybutbreaksthetime-reversalsymmetry. Mean- − x origin k = 0 in the wavevector space within the WSM while,wecanalsodescribeWSMswithoutinversionsym- phase. Finally, in the limit of a/b the Weyl nodes metry by combining two models described by Eq. (1), → ∞ which is beyond the scope of the present paper. converge to (±√m+,0,0). Thus the multilayer is in the WSM phase when the NI layer is sufficiently thin com- The two Weyl nodes in this model are displaced along paredto the WSM layer. Fromthese considerations,the the kx direction in the wavevector space (Fig. 1 (a)). phase diagram is as shown in Fig. 2 (b). Therefore, in multilayers of a WSM and a NI, the stacking direction relative to the displacement vector w between the two Weyl nodes (i.e. k direction) will affect x their properties. Here, we study two cases for the multilayers: Patterns A and B. In Pattern A (Fig. 1 (b)), the stackingdirectionisperpendiculartothevectorw,while in Pattern B (Fig. 1 (c)), it is parallel to the vector w. B. Multilayer - Pattern A In this subsection, we discuss the WSM-NI multilayer of Pattern A using the effective model (1). The stacking direction is in the z direction, perpendicular to the dis- FIG. 1. (a) Schematic diagram of the Weyl nodes, displaced placement vector w between the two Weyl nodes. Since along the kx direction in the wavevector space. The vector the model (1) represents a WSM when m > 0 and a NI w is the displacement vector between the two Weyl nodes in when m < 0, the multilayer is realized by putting the thekspace. (b),(c)Stackingpatternsofthemultilayers: (b) Pattern A and (c) Pattern B. control parameter m to be periodic as shown in Fig. 2 3 We calculate the eigenstates and the energies in the similar way as in the Krönig-Penney model. The details of the calculation are given in Appendix A2. In order to obtain the phase diagram as a function of a and b, we examine whether the band gap closes. The band gap closes when m m − + − sin√m+asinh√m−b+cos√m+acosh√m−b 2√m+m− =cosk (a+b), k =k =0 (13) x y z FIG. 2. (a) Schematic diagram of the multilayer of Pattern is satisfied. Here, k is the Bloch wavenumber along the x A and the z dependence of the parameter m(z). (b) Phase x axis. Solutions of Eq. (13) give positions of the Weyl diagram of the WSM-NImultilayer of Pattern A. nodes ( k ,0,0). When Eq. (13) has (real) solutions of x ± k , the system is in the WSM phase. This WSM phase x is in the orange regions in Fig. 3 (b) for m = 0.5 and Next, we calculate the surface state of this multilayer. + m =0.25 as an example. − We consider the multilayer with semi-infinite geometry In the other regions, the bulk is gapped. Because the with the z =0 surface;the z <0 regionis the multilayer time-reversalsymmetryis broken,these regionswith the andthez >0regionisthevacuum. Thewavefunctionof bulk gap can be either the NI phase or the QAH phase. the surfacestateanditsenergyeigenvaluearecalculated One can determine the phases in these regions in the as followingway. Bycontinuity,thephasesintheindividual ψ = 1 e−(γ/v)Rzdz(kx2−m(z)), E = vk . (7) regionsinthephasediagramcanbeeasilydeterminedby i − y consideringthelimitb withfixeda,asshowninthe (cid:18)− (cid:19) →∞ Thisdescribesthesurfacestateonlywhentheintegralin Eq. (7) decays into the bulk; this condition is given by (a+b)k2 (m a m b)<0. (8) x− + − − Hence, the surface state exists when m a m b m a m b + − + − − <k < − ; (9) x − a+b a+b r r namely,itappearsbetweenthetwoWeylnodesandforms a Fermi arc. In addition, the surface state has a velocity (0, v) along the y direction from Eq. (7). From the − bulk-edge correspondence, we conclude that the Chern number of the two Weyl nodes m+a−m−b,0,0 is ± a+b 1, respectively. (cid:18) q (cid:19) ± C. Multilayer - Pattern B In this section, we consider the WSM-NI multilayer of PatternB.Thestackingdirectionisalongthexdirection, FIG. 3. (Color online) (a) Schematic diagram of the multi- which is parallel to the displacement vector w between layer of Pattern B and the x dependence of the parameter the Weyl nodes. We set the parameter m to be periodic m(x). (b) Phase diagram by changing the thicknesses a and as shown in Fig. 3 (a), and the Hamiltonian is given by b in the multilayer when m+ = 0.5, m− = 0.25. The WSM phaseisintheorangeregionsandtheNIphaseisinthewhite H =γ −∂x2−m(x) σx+v(kyσy+kzσz), (10) region. EachblueregionistheQAHphasewithdifferentval- ues of the Chern number ν. (c) Positions of the Weyl nodes (cid:0) m+ (0(cid:1) x a), m(x)= ≤ ≤ (11) asthethicknessaischanged,withb=1fixedandm+ =0.5, ( −m− (−b≤x≤0), m− =0.25. TheWeylnodesmovealongthekx axis. Thered linerepresentsthetrajectory ofamonopoleandtheblueline m(x+(a+b))=m(x), (12) representsthat of an anti-monopole. They correspond tothe change of a along thedashed line in (b). where m are positive constants. ± 4 Appendix B. In this limit, the system reduces to a thin slab of a WSM, and one can easily calculate its Chern number. We then obtain the phase diagram Fig. 3 (b). Here, the QAH phases are characterized by the Chern number ν within the constant k plane. The values of z the Chern number of two phases separatedby the WSM phase are different by one. WealsoinvestigatemovementoftheWeylnodesinthe wavevectorspacewhenthethicknessoftheWSMlayeris changed. Here,wefixthethicknessoftheNIlayerasb= 1, and we gradually increase the thickness of the WSM layera. ThepositionsoftheWeylnodesalongthek axis x areshowninFig.3(c)whenm =0.5,m =0.25. First, + − whenthemultilayerenterstheWSMphaseasaincreases FIG. 4. (Color online) (a) Three-dimensional unit cell of the fromzero,apairofWeylnodesiscreatedatk=0. Then, lattice model. (b)Two-dimensional unit cell in they-z plane the Weyl nodes continuously move along the k axis as (c) Part of the unit cell along the (110) and (1¯10) planes, x a increases, and finally they are annihilated pairwise at whichcontaintheAandBsub-lattices. Thesolidandbroken the boundary of the Brillouin zone in the x direction. lines represent the nearest-neighbor hopping and the next- Such pair creations and annihilations alternately occur nearest-neighbor hopping, respectively. Colored regions in as we increase a further. Therefore, the phases with a (b)and(c)arethreadedbyamagneticflux: φ0/4intheblue bulk gap and the WSM phase alternately appear in the regions,andφ0/2intheorangeregion,whereφ0 =h/eisthe magnetic flux quantum. phase diagram of the multilayer. From the increment of the value of the Chern number by the increase of a, one can identify which trajectory is of a monopole or of an The Hamiltonian of the lattice model is written as anti-monopole; the result is shown in Fig. 3 (c). H =(2tsink a )σ +(2tsink a )σ + x − y · · +[∆ 2t′(cosk a +cosk a )+2t cosk a ]σ , y z ⊥ x z − · · · (14) III. MULTILAYER FROM THE LATTICE MODEL where a , a and a are primitive lattice vectors, and x y z a =(a +a )/2. Here, we choose the coordinate axes ± y z In this section, we use the lattice model for a WSM asshowninFig.4(a), whicharedifferentfromRef.[32], proposed in Ref. 32 in order to study the WSM-NI mul- in order to set the displacement vector w between the tilayers. InSec.IIIA,weintroducethelatticemodelfora Weyl nodes to be along the x axis in accordance with WSM.Byusingthislatticemodel,wenumericallycalcu- Fig. 1 (a). For simplicity, we set ax = ay = az = d. | | | | | | latethebandstructureandphasediagramsforthemulti- The bulk energy is given by layerswithPatternsAandBinSec.IIIBandSec.IIIC, E = 2 t2 sin2k d+sin2k d respectively. ± + − 1/2 + t2 ((cid:2)m(cid:0) m cosk dcosk (cid:1)d+cosk d)2 (,15) ⊥ 1− 2 + − x i where k =(k +k )/2, m =∆/t and m =2t′/t . ± y z 1 ⊥ 2 ⊥ A. WSM from the lattice model The bulk gap closes when sink d = sink d = 0 and + − cosk d = (m m cosk dcosk d). Within the re- x 1 2 + − − − In this subsection, we review the lattice model pro- giongivenbym +m >1, m m <1andm ,m > 1 2 1 2 1 2 posed in Ref. 32. It is a tight-binding model of spinless 0, the gap closes at W =| −1arcc|os(m m ),0,0 ± ±d 2− 1 fermions on a lattice of stacked face-centered squares, in the three-dimensional wavevector space. This region (cid:0) (cid:1) shown in Fig. 4 (a). The unit cell of the lattice model represents the WSM phase. On the other hand, when includessub-latticesAandB.Thenearest-neighborhop- m m > 1 and m , m > 0, the bulk is gapped, and 1 2 1 2 − pingbetweentheAandBsub-latticeshasanamplitudet, the system is the NI phase. andthenext-nearest-neighborhoppingbetweenthesame In the following subsections Sec. IIIB and Sec. IIIC, sub-lattices has an amplitude t′. Furthermore, the hop- we numerically calculate the bulk bands of the WSM-NI pingamplitudealongthestackingdirectionisdenotedby multilayersby use ofthe lattice model(14)andcompare t between the same sub-lattices on adjacent layers. We the results with those in Sec. IIB and Sec. IIC. In order ⊥ also add an on-site energy for the A and B sub-lattices: torealizetheWSM-NImultilayerusingthelatticemodel +∆ for the A sub-lattice and ∆ for the B sub-lattice. (14), the parameters m , m are periodically modulated 1 2 − The magnetic flux is added to the lattice with the gauge between those for a NI and those for a WSM. For sim- choiceshowninFigs.4(b)and(c)withφ =h/e,where plicity, we fix the value of m and change the value of 0 2 h is the Planck constant and e is the electron charge. m to be m in the WSM layer and m in the NI 1 1W 1N − 5 layer, where the parameters m and m should meet adjacent Fermi arcs through the WSM layer and those 1W 1N the above conditions. Furthermore, we set N and N through the NI layer,respectively. Furthermore, let L a b ≡ to be the numbers of the atomic layers within the WSM a+b denote the period of the multilayer. An effective layer or the NI layer, respectively. Then, the thickness modeldescribinghoppingbetweentheFermiarcsisthen of the WSM layer and that of the NI layer are given by written as a=Nad and b=Nbd, respectively. 0 t+t′e−ikzL H = , (16) t+t′eikzL 0 (cid:18) (cid:19) B. Multilayer - Pattern A and its energy eigenvalues are given by E = t2+t′2+2tt′cosk L. (17) ± z In this subsection, we study the multilayer of the lat- ± tice model stacked along the z direction, which is per- Thehoppingamplitupdest,t′dependonk andk . Phys- x y pendicular to the directionof the displacement vector w ically, they depend on magnitudes of the band gaps of between the two Weyl nodes. the two layers. For example, by increasing the thickness We numerically calculate the phase diagram of the of the WSM layer or by increasing the band gap of the multilayerfromthelatticemodelwiththeparameterval- WSM layer, the hopping amplitude t through this layer ues m = 1.5, m = 2.5, m = 1.0 as shown in Fig. 5 asymptoticallybecomezero. Thesimilarbehaviorisseen 1w 1N 2 (a). In the phase diagram, the orange circles represent also for the NI layer. Namely, for the multilayer, t′ 0 → the WSM phase and the green squares represent the NI when we increase the thickness or the band gap of the phase. This phase diagram qualitatively agrees with the NI layer. On the other hand, t 0 when we increase → result of Fig. 2 (a) from the effective model. Examples thickness of the WSM layer at a wavenumber where the of the band structure are shown in Fig. 5 for N = 5, WSM layer has a gap. In order to confirm this scenario, a N =5 in (b) and for N =100, N =5 in (c). In these we investigatemagnitude ofthe hopping amplitudes t, t′ b a b cases,theWeylnodes,W andW ,appearonthek axis as functions of the numbers of the layers in the multi- 1 2 x (k = k = 0). Comparing (b) and (c), we find that the layer using Eq. (17). For the hopping amplitude t with y z positionsoftheWeylnodesmoveinthewavevectorspace the parameter values m = 1.5, m = 2.1, m = 1.0, 1W 1N 2 by changing the thickness of the WSM layer a = N d. N = 30, k = 0, k = π, the result shown in Fig. 6 (c) a a y x It corresponds to the result in Eq. (6). Namely, as we fitswellwithexponentiallydecayingformt exp( n/λ) decrease N in Fig. 5 (a), the Weyl nodes gradually ap- with λ = 5.83. Similarly, the hopping amp∝litude t−′ with a proach k = π. Further decrease of N causes a pair m = 1.5, m = 2.5, m = 1.0, N = 30, k = 0, x a 1W 1N 2 b y annihilationo±fWeylnodesatk =π,drivingthe system k = 2.1 fits well with t′ exp( m/λ′) with λ′ = 10.0. x x ∝ − into the NI phase. These Weyl nodes affect the Chern Thus,this picture wellexplainsthe statesclose to E =0 numberonthek k plane. ThedependenceoftheChern in Fig.6 (b). According to Eq.(17), the bandgapcloses y z number ofthe wavenumberk is showninFig.5(d) and when t(k ,k ) = t′(k ,k ) within the WSM phase, x x y x y ± (e) for N = 5, N = 5 and N = 5, N = 100, re- and it is realized at some k with k =0. Meanwhile, it a b a b x y spectively. The Chern number on the plane k = const. is not satisfied in the NI phase and there appears a gap x changesby 1attheWeylnodes. Fromthesefigures,we 2 t t′ . ± | ± | conclude that the Weyl node W has a monopole charge 1 +1, and the Weyl node W has a monopole charge 1. 2 − Intheabovediscussion,wefoundthatthemultilayeris C. Multilayer - Pattern B intheWSMphasewhenthe WSMlayeris thickenough. Ontheotherhand,whentheNIlayeristhick,itisinthe We calculate the band structure for the multilayer of NI phase. Nevertheless, for a sufficiently thick NI layer, PatternB.Weshowournumericalresultsofthephasedi- the band structure, shown in Fig. 6 (b), looks quite dif- agramand the Chern number with changing thicknesses ferent from that for the NI phase with uniform value of ofthe WSM andthe NI layersinFig.7. Inthe phasedi- m = m as shown in Fig. 6 (a). Namely, it is simi- agram,the regionsshownby the orangecirclesrepresent 1 1N lar to the band structure of the bulk of the NI (Fig. 6 the WSM phase, and the regions shown by the squares (a)),butwithadditionalbandsaroundtheFermienergy. representphaseswithabulkgap. Here,atthepointsrep- This behavior is caused by the Fermi arc which appears resented by the white circles, we could not identify the on the interfaces between the WSM layers and the NI phasesnumerically. Thephaseswithabulkgap(squares layers. Namely, since the adjacent Fermi arcs are sepa- in Fig. 7) are further classified as follows. Green repre- rated by the insulating layer, the states near the Fermi sentstheNIphase,andlightblue,purpleandbluerepre- energy come from the Fermi arcs on the interfaces, with sentthe QAHphasewiththe Chernnumber 1, 2and − − small hybridization. In Fig. 6 (b), there is a tiny gap 3,respectively. TheChernnumberonthek k planeas y z − throughout the whole Brillouin zone, and it represents a function of k for variousvalues of the thickness of the x a NI phase. In the following, we show this by a simple WSM layer is shown in Fig. 7 for (b) N = 3, N = 2, a b model. (c) N = 4, N = 2 and (d) N = 5, N = 2. Fig- a b a b Lettandt′denotethehoppingamplitudesbetweenthe ure. 7 (b) (N =3,N =2) represents the QAH phase a b 6 FIG. 5. (Color online) (a) Phase diagram of the multilayer from the lattice model with the parameter values m1W = 1.5, m1N = 2.5, m2 = 1.0. The orange circles represent the WSM phase and the green squares represent the NI phase. (b), (c) Bulk band structure and the Chern number of the multilayer with the parameter values m1W = 1.5, m1N = 2.5, m2 = 1.0. The thicknesses of the layers are (b) Na = 5, Nb = 5 and (c) Na = 100, Nb = 5. The Weyl nodes move along the kx axis by changing the thickness of the WSM layer a = Nad. (d), (e) Chern number of the multilayer with the parameter values m1W =1.5, m1N =2.5, m2=1.0 and (d) Na =5, Nb =5 and (e) Na =100, Nb =5. with ν = 1 for any value of k . As the thickness of ically appears as the thickness of the WSM layer varies, x − the WSM layer becomes larger as shown in Fig. 7 (c) and two phases which sandwich the WSM phase are the (N =4,N =2), the multilayer is in the WSM phase. QAH phases with different values of the Chern number. a b Here, the Weyl nodes divide the k axis into the regions Agreement between the results from the effective model x of ν = 1 and of ν = 2. The phase transition between and those from the lattice model is surprisingly good. − − Fig.7(b)and(c)isattributedtoapaircreationofWeyl The band structure in Pattern A with relatively thick nodes,andweinterpretthatitoccursattheboundaryof layers (e.g. Fig.6 (b)) can be understood from mixing of the first Brillouin zone. As Na becomes larger,the Weyl theFermiarcsattheindividualinterfaces,withsmallhy- nodes move further, until they annihilate pairwise. The bridization between them. The exponential dependence multilayer is then in the QAH phase with ν = 2 for ofthehybridizationasafunctionofthethickness(Fig.6 − allkx (Fig.7(d)). Therefore,twobulk-insulatingphases (c)) confirms this scenario. Furthermore, we can apply whichsandwichthe WSMphase havethe Chernnumber the scenario of the mixing of the Fermi arcs not only in different by 1. This conclusion agrees well with that PatternAbutalsoinPatternB.InPatternB,theFermi − fromthe effective model. Therefore,we haveshownthat arcs in the WSM layer are folded down into the first the WSM-NI multilayer undergoes the phase transitions Brillouin zone due to the periodic stacking. It gives rise with various phases including the WSM phase and the to multiple chiral surface states circulating along the yz QAH phase with various Chern numbers. plane. As a result, the QAH phase appearson the phase diagrams (Fig. 3 (b) and Fig. 7 (a)). Therefore, appearance of the QAH phase can be regarded as a result of IV. SUMMARY AND DISCUSSION the hybridization between the Fermi arcs. Interestingly, the Chernnumber ofthe QAHphasescanbe largewhen Inthispaper,westudybehaviorsoftheWSM-NImul- we increase the thickness of the WSM layer a, and with tilayerusingtheeffectivemodelandthelatticemodel. In this multilayer, one can design a QAH phase with large both models, the behaviors of the multilayer are differ- Chern number, which has been a long-standing issue in ent for different stacking directions. In the multilayer of this field33–38. Pattern A, i.e. with the stacking direction perpendicu- To the authors’ knowledge, WSM-NI multilayers have lar to the displacement vector w of the Weyl nodes, the not been studied previously. On the other hand, the phase diagram is obtained as Fig. 2 (a), and the multi- WSMphaseisfoundtoappearinaTI-NImultilayerwith layer undergoes the phase transition from the NI phase magnetization26andinastakedQAHlayerssystem39. In to the WSM phase by increasing the thickness of the the phase diagrams of these two cases (TI-NI multilayer WSM layer. On the other hand, in PatternB,where the and stacked QAH layers system), the WSM phase ap- stacking direction is parallel to the displacement vector pearsbetweenthe NIandthe QAHphases. Bychanging w of the Weyl nodes, the phase diagram is obtained as system parameters, one can go through NI-WSM-QAH Fig.3(a). Inthephasediagram,theWSMphaseperiod- phase transitions. Here, the WSM can be seen as an in- 7 FIG. 7. (Color online) (a) Phase diagram for the multilayer for the parameter values m1W = 1.5, m1N = 2.1, m2 = 1.0. The circles represent the WSM phase, and the squares represent the phases with a bulk gap. The colors of the symbol referas differentphases. (b)-(d)Chern numberas afunction of kx with (b) Na = 3, Nb = 2, (c) Na = 4, Nb = 2 and (d) Na = 5, Nb = 2. (b) shows the QAH phase with the Chern FIG. 6. (a) Bulk band structure of the lattice model in the number −1 for all kx, (c) shows the WSM phase, and (d) NI phase for the parameter values Nb = 100, m1N = 2.5, shows the QAHphase with theChern number−2 for all kx. m2 =1.0. (b)Band structurefortheWSM-NImultilayerfor the parameter values Na =5, Nb =100, m1W =1.5, m1N = 2.5, m2 = 1.0. A band gap is almost zero as shown as red vectorn, defined asa normaldirectionofindividuallay- lines. (c) Hopping amplitudes across the layers as a function ers,iseitherperpendicularorparalleltothedisplacement of thenumbersof atomic layers, shown in asemi-logarithmic′ vector w between the two Weyl nodes. As we have dis- plot. The red squares represent the hopping amplitude t cussed, in Pattern A (n w), the Weyl nodes approach through the NI layer with the parameter values m1W = 1.5, ⊥ each other by making a multilayer, and in Pattern B m1N = 2.1, m2 = 1.0, Na = 30, ky = 0, kx = π. The (n w), the k space is folded into the Brillouin zone by green circles represent the hopping amplitude t through the k WSMlayerwiththeparametervaluesm1W =1.5,m1N =2.5, the periodicity of the multilayer and becomes the WSM, m2 =1.0, Nb =30, ky =0, kx =2.1. The red and the green the NI or the QAH phases depending on the thicknesses lines show a fittingby exponentialfunction. of the two layers. In the WSM with N = 2, one can think of general multilayers, where the two vectors n and w are neither complete version of the QAH; namely, at some 2D slices parallel nor perpendicular. Such general cases can be in k space the Chern number is zero,while at others the understood by a combination of Patterns A and B. In Chern number is nontrivial. In this sense, our finding of the directionperpendicular to the stackingvectorn, the an appearance of the QAH phase in the WSM-NI multi- Weylnodesapproacheachother,whilealongthestacking layer is nontrivial. As stated in the previous paragraph, vector n, the k space is folded into the small Brillouin theemergentQAHphaseinPatternBisunderstoodfrom zone. Thus the displacement vector between the Weyl thefoldingoftheBrillouinzonebyperiodicityofthemul- nodes becomes gradually parallel as we increase b/a. As tilayer. long as n is not parallel to the displacement vector be- In this paper, we consider the WSM with broken tweentheWeylnodes,the Weylnodeshavenochanceof time-reversal symmetry. There are a number of ma- pair annihilation, and the system is in the WSM phase. terials proposed as WSMs with broken time-reversal Then, when b/a exceeds a critical value, as is similar to symmetry: pyrochlore iridates5,8,9, CoTi X (X=Si, Ge, Pattern A, the displacement vector between the Weyl 2 Sn)13 ZrCo Sn14, Hg Cd Mn Te in a magnetic nodes becomes parallelto n, and the multilayermay ex- 2 1−x−y x y field18, and Bi Sb in a magnetic field19. The hibittheQAH,theWSM,ortheNIphases. Thedetailed 0.97 0.03 number of Weyl nodes, N, takes various even num- behavior of the phases as a function of the thicknesses a bers in these WSMs. In WSMs with N = 2, and b depends on details of the system, and is left as a such as Hg Cd Mn Te under a magnetic field18, future work. 1−x−y x y HgCr Se 10,11,andBi Sb underamagneticfield19, In the other WSMs with brokentime-reversalsymme- 2 4 0.97 0.03 thereareonlyapairofWeylnodesaroundthe Fermien- try having N >2, the phase diagram for the multilayers ergy. The presenttheoryisdirectly applied,ifthe multi- canbediscussedinasimilarway,andtheresultingphase layer is either of the two patterns A and B; i.e. stacking diagram will become complicated. Nevertheless, based 8 on the above analysis, one can say that possible phases where the A, B, C, D, K and Q are constants. are either a NI, a WSM, or a QAH phases, because the Hence, by continuity of the wave function and the layered structure folds the Fermi arcs into the Brillouin Bloch condition, we obtain zone of the multilayer. In general, a longer periodicity in the superlattice direction gives rise to larger Chern ψ(+0)=ψ( 0), − (A3) number, because the Fermi arcs are folded down many ψ(a 0)=ψ( b+0) eikz(a+b), (cid:26) − − times. An extension to the WSM without the inversion sym- wherekz istheBlochwavenumberalongthez direction. metry is beyond the scope of the present work, as the By combining Eqs. (A1) - (A3), we obtain the condition symmetry of the systems is very different. There are for the coefficients A, B, C and D to have non-trivial at least four Weyl nodes, and when the folding of the values: Brillouin zone gives rise to pair annihilations of all the 1 Q2 K2 γ2(m +m )2 Weyl nodes, it is expected to lead to either a NI or a − + − sinKasinhQb weak TI phases. Nevertheless, the way how the Fermi KQ" 2 − 2v2 # arcs connect the four Weyl nodes depends on surface +cosKacoshQb=cosk (a+b). (A4) termination31; therefore, a discussion from a viewpoint z of Fermi arcs is not straightforward, and is left as a fu- SinceKandQarefunctionsoftheenergyEviaEq.(A2), ture work. Eq. (A4) determines the energy eigenvalues as functions of the Bloch wavevector k=(k ,k ,k ). x y z Next, we derive the condition for closing of the band ACKNOWLEDGMENTS gap. We have set the Fermi energy to be E = 0. Be- F cause the energy bands are symmetric with respect to ThisworksupportedbyGrant-in-AidforScientificRe- E = 0, the gap closing condition is obtained by setting search(No. 26287062,16K13834)byMEXT,Japan,and E = 0 in Eq. (A4). This yields Eq. (6) after a straight- byMEXTElementsStrategyInitiativetoFormCoreRe- forward calculation. search Center (TIES). 2. Multilayer - Pattern B Appendix A: Calculations of eigenstates and energies for the WSM-NI multilayers Inthissubsection,weexplainthedetailsofcalculations of eigenstates and energy eigenvalues for the WSM-NI 1. Multilayer - Pattern A multilayer with Pattern B in Sec. IIC. First, we write downthe wavefunctionandthe energyeigenvalueinthe Here, we explain the details of calculations of eigen- WSM and the NI layers. In the WSM layer, the wave statesandenergyeigenvaluesfortheWSM-NImultilayer function is given by withPatternAinSec.IIB.First,wewritedownthewave function and the energy eigenvalue both in the WSM γ K2+m +ivk layer and in the NI layer, and connect the wave func- ψ(x)= − 1 + y A1eiK1x+B1e−iK1x vk E tion at the boundaries. The wave function ψ(z) in the (cid:18) (cid:0) z − (cid:1) (cid:19) WSMlayer(0<z <a)andintheNI layer(−b<z <0) + γ −K22+m+ +ivky (cid:0)C1eiK2x+D1e−iK2x(cid:1) is given respectively by vk E (cid:18) (cid:0) z − (cid:1) (cid:19) (cid:0) (cid:1) E+vK (0<x<a), (A5) ψ(z)=A eiKz γ k2 m +ivk (cid:18) x− + y (cid:19) and the energy eigenvalue is written as E vK +B (cid:0) − (cid:1) e−iKz (0<z <a), γ k2 m +ivk E = γ2(K2 m )2+v2 k2+k2 (cid:18) x− + y (cid:19) ± 1 − + y z E ivQ q ψ(z)=C (cid:0) − (cid:1) eQz = γ2(K2 m )2+v2(cid:0)k2+k2(cid:1), (A6) γ k2+m +ivk ± 2 − + y z (cid:18) x − y (cid:19) q (cid:0) (cid:1) +D (cid:0) E+iv(cid:1)Q e−Qz ( b<z <0), where A1, B1, C1, D1, K1 and K2 are constants. The γ k2+m +ivk − wave function and the energy eigenvalue in the NI layer (cid:18) x − y (cid:19) ( b<x<0) are given by replacing m m , K + − i (A1) − →− → (cid:0) (cid:1) iQ (i=1,2), A A , B B , C C and i 1 2 1 2 1 2 − → → → and the energy eigenvalue E is written as D1 D2. Here, we select K1 and Q1 as → E =± γ2(kx2−m+)2+v2 ky2+K2 K12 =m++ γ1 E2−v2 ky2+kz2 , (A7) =±qγ2(kx2+m−)2+v2(cid:0)ky2−Q2(cid:1), (A2) −Q21 =−m−+qγ1 E2−(cid:0)v2 ky2+(cid:1)kz2 q q (cid:0) (cid:1) (cid:0) (cid:1)  9 thickness a, and this slab can be described within the effective model, when the parameter m is as shown in Fig. 8. Let us calculate the Chern number for this WSM slab with the thickness a in Fig. 8. Here, the Chern number refers to that defined in the k k plane. By the mir- y z ror symmetry with respect to the yz plane, the mirror eigenvalues of the wave functions are either = +1 or M = 1, yielding ψ( x) = ψ(x) and ψ( x) = ψ(x), M − − − − respectively. Namely, the wavefunctions are either symmetric or antisymmetric with respect to x=0. FIG.8. Valueofm(x)fortheslaboftheWSMwiththickness Next,wewritedownthewavefunctionsineachregion. a. It corresponds to themultilayer in the limit b→∞. The wavefunctions in the regiona/2<x canbe written as and K2 and Q2 as ψ(x)= γ Q21−m− +ivky Ae−Q1x vk E · K2 =m 1 E2 v2 k2+k2 , (cid:18) (cid:0) z −(cid:1) (cid:19) 2 +− γ − y z (A8) + γ Q22−m− +ivky Be−Q2x a <x ,  Q2 = m q1 E2 (cid:0)v2 k2+(cid:1)k2 . vk E · 2  − 2 − −− γ − y z (cid:18) (cid:0) z −(cid:1) (cid:19) (cid:16) (cid:17) q (B1) By continuity ofthe wavefunction a(cid:0)ndthe B(cid:1)lochcon- dition we obtain and the energy eigenvalues are given by ψψ((a+0)0=)=ψ(ψ−(0)b, +0) eikx(a+b), E =± γ2(−Q21+m−)2+v2 ky2+kz2  − − (A9) q  ψψ′′((a+0)0=)=ψ′ψ(−′(0)b,+0) eikx(a+b), =±qγ2(−Q22+m−)2+v2(cid:0)ky2+kz2(cid:1). (B2) − − Here A, B, Q and Q are constants(cid:0). Next,(cid:1)we write wherek istheBlochwavenumberalongthexdirection. down the wave1function2 in the region a/2 < x < a/2. x − Then, the condition for the coefficients A , B , C and The wave functions with a mirror eigenvalue = +1 i i i M D (i=1,2) to have non-trivial values is are given by i Q22iK−QKi2 sinKiasinhQib+cosKiacoshQib ψ(x)= γ −K12v+kzm+E+ivky ·CcosK1x i i (cid:18) (cid:0) − (cid:1) (cid:19) =coskx(a+b) (i=1,2). (A10) + γ −K22v+kzm+E+ivky ·DcosK2x Since K andQ (i=1,2)arefunctions ofthe energyE, (cid:18) (cid:0) − (cid:1) (cid:19) i i a a Eq.(A10)determinestheenergyeigenvaluesasfunctions <x< , (B3) of the Bloch wavevectork=(k ,k ,k ). −2 2 x y z (cid:16) (cid:17) We can derive the condition for closing of the band andthosewithamirroreigenvalue = 1aregivenby M − gap by setting E = 0 in Eq. (A10), because the energy γ K2+m +ivk bands are symmetric with respect to E =0. This yields ψ(x)= − 1vkz +E y ·CsinK1x Eq. (13) after a straightforwardcalculation. (cid:18) (cid:0) − (cid:1) (cid:19) γ K2+m +ivk + − 2vkz +E y ·DsinK2x (cid:18) (cid:0) − (cid:1) (cid:19) Appendix B: Calculation of the Chern number of a a a <x< . (B4) WSM slab from the effective model −2 2 (cid:16) (cid:17) Here, the energy eigenvalue is given by AsseenfromFig.3(a),thephasediagramforPattern B in the region with a large value of b consists mostly E = γ2(K2 m )2+v2 k2+k2 of bulk-insulating phases, separated from each other by ± 1 − + y z q the WSM phase, having a narrow region. These bulk- = γ2(K2 m )2+v2(cid:0)k2+k2(cid:1), (B5) insulating phases mightbe the NI with zeroChern num- ± 2 − + y z q ber or the QAH phase with nonzero Chern number. To and C, D, K1 and K2 are constants,(cid:0)and we(cid:1)select K1, calculatethe Chernnumberforeachphase,wenote that Q1 as Eq. (A7) and K2, Q2 as Eq. (A8). the systeminthisregionofthe largevalueofbisasymp- Now, by the continuity condition, totically a collection of WSM layers, separated far away ψ a +0 =ψ a 0 , from each other. Thus, the Chern number can be cal- 2 2 − (B6) culated by considering a single slab of a WSM with the (ψ′(cid:0) a2 +0(cid:1) =ψ(cid:0)′ a2 −(cid:1)0 , (cid:0) (cid:1) (cid:0) (cid:1) 10 we derive a set of equations for the coefficients A, B, C, D, which eventually decouples to equations for A and C, and those for B and D. In either of these coupled equations,wecanproceedinthesimilarwaybychanging the notation K K, Q Q or K K, Q Q. 1 1 2 2 → → → → Here,weadoptB =D =0withoutlosinggeneralityThe energy eigenvalues of the system are given by E2 =γ2M2+v2 k2+k2 , (B7) y z (cid:0) (cid:1) where M is defined from K, Q as M K2 m = Q2+m . (B8) + − ≡ − − Finally, when the mirror eigenvalue is = +1, the M condition for the existence of non-trivial solutions from FIG. 9. (Color online) Dependence of the parameter Mn on Eq. (B6) is given by the slab thickness a determined by Eq. (B9) (red lines) and Eq.(B10)(bluelines)whenm+ =0.5andm− =−0.25. The Chern numberν of the system changes by unityat thevalue m +Ma m M tan + = −− . (B9) of a such that Mn =0. p 2 sm++M where When a mirror eigenvalue is = 1, the condition is given by M − ∂ψ− ∂ψ− ∂ψ− ∂ψ− nk nk nk nk ∂k ∂k − ∂k ∂k (cid:28) y (cid:12) z (cid:29) (cid:28) z (cid:12) y (cid:29) tan pm+2+Ma − π2!=smm−+−+MM. (B10) =Z−∞∞dx(cid:12)(cid:12)(cid:12) ∂ψ∂n−kk†y(x)∂ψ∂n−kkz(x)(cid:12)(cid:12)(cid:12)− ∂ψ∂n−kk†z(x)∂ψ∂n−kky(x)!, (B14) Thus, from Eqs. (B9) and (B10), the thickness of the WSM layer a determines the values of the parameter and ψ− is the nth eigenstate at the Bloch wavevector nk M =Mn as k = (ky,kz). The summation of the right-hand side of Eq. (B13) is taken over the states ψ− below the Fermi nk m +M a nπ m M energy. Here, the integral is over the two-dimensional + n − n = +arctan − , p 2 2 sm++Mn fBercitlilvoeuimnozdoenlediensctrhibeeksytkhzepsylasnteem. Wonelynnoetearthka=t t0h,eanefd- n=0,1,2, . (B11) ··· we assume that the system itself has the periodicity of the Brillouin zone. As we increase the thickness a, the For each value of n, the energy eigenvalues are given by Chern number changes only when the band gap closes, Eq. (B7): i.e. M = 0 for some n; this occurs at the intersections n betweenthe curvesandthe horizontalaxisinFig.9. Let En± ≡± γ2Mn2+v2 ky2+kz2 . (B12) aTh=eajunemdpenooftethteheChthericnknneusmsbwehrearet aM=ne =a 0isfocralscoumlaetend. q n (cid:0) (cid:1) as e Thus,nand serveasabandindex,andthecorrespond- e ingeigenstat±esψ± haveamirroreigenvalue =( 1)n. 1 n M − ν(a +0) ν(a 0)= ∆sgn(M ) (B15) Inparticular,fromEq.(B12),thesystemisgaplesswhen n n n M = 0. We show the dependence of M on the slab e − e− 2 e (cid:12)a=an n n (cid:12) e thickness a in Fig. 9. At the thickness a such that where the right hand side is 1 when Mn(cid:12)(cid:12)changes from Mn =0, the gap of the nth bands vanishes. positivetonegativeata=an−,andis+1wehenitchanges Next, we calculate the Chern number of the system from negative to positive asewe increase a . Therefore, n from the wave functions calculated above. The Chern from Fig. 9, the Chern number changes bye 1 as a in- − number on the k k plane is defined as creases across a . At a = 0, the Chern number of the y z n system is zero sience the system is in the insulator phase. ∞ dkydkz As a gradually increases, the parameter Mn=0 becomes ν = zero. TheChernnumberofthesystemthenchangesfrom 2π nX=0ZB.Z. zero to 1. As a increases further, Mn=1 becomes zero i ∂ψn−k ∂ψn−k ∂ψn−k ∂ψn−k ,(B13) andthe−Chernnumberofthesystemchangesfrom−1to × ∂k ∂k − ∂k ∂k 2. Thus, we have determined the Chernnumber for all (cid:18)(cid:28) y (cid:12) z (cid:29) (cid:28) z (cid:12) y (cid:29)(cid:19) − (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12)