Electrically tunable multiple Dirac cones in thin films of (LaO)2(SbSe2)2 family of materials

Two-dimensional Dirac physics has aroused great interests in condensed matter physics ever since the discovery of graphene and topological insulators due to its importance in both fundamental physics and device applications. The ability to control the properties of Dirac cones, such as bandgap and Fermi velocity, is essential for the occurrence of various new phenomena and the development of next-generation electronic devices. Based on first-principles calculations and an analytical effective model, we propose a new Dirac system with eight Dirac cones in thin films of the (LaO)2(SbSe2)2 family of materials with an external gate voltage. The advantage of this system lies in its tunability: the existence of gapless Dirac cones, their positions, Fermi velocities and anisotropy all can be controlled by an experimentally feasible gate voltage. We identify the layer dependent spin texture induced by spin-orbit coupling as the underlying physical reason for the tunability of Dirac cones in this system. As a consequence, we show that the electrically tunable quantum anomalous Hall effect with a high Chern number can be induced by introducing magnetization into this system.


I. INTRODUCTION
Unlike the usual materials with parabolic energy dispersion described by the non-relativistic Schrödinger equation, graphene 1 provides the first example in condensed matter physics with low energy effective physics described by the relativistic Dirac equation with linear energy dispersion. Later, it was realized that two dimensional Dirac type of dispersion also appears for the surface states of three dimensional topological insulators (TIs) 2,3 , of which spin is resolved and locked to the momentum, forming a spin texture in the surface Brillouin zone (BZ). Dirac Hamiltonians with or without mass also exist in the low energy physics of several other families of materials, including topological crystalline insulators (such as SnTe 4-6 ), group-VI dichalcogenides (such as MoS 2 7-11 ), SrMnBi 2 12 , d-wave cuprate superconductors 13 , and three-dimensional Dirac semimetals (such as Na 3 Bi 14 , Cd 3 As 2 15,16 ). Due to the similar energy dispersion, all these materials share some common and unique physical properties, thus dubbed "Dirac materials" 17 , which are believed to have potential applications in high performance nanoelectronics 1 , spintronics, and quantum computation 2,3 . In these existing Dirac materials, the properties of Dirac cones, such as the position and Fermi velocity, are usually determined by intrinsic properties of material, such as crystal structures and spin-orbit coupling (SOC), and thus difficult to be controlled experimentally. For example, pristine graphene 1 does not have a bandgap, so for the potential application in transistors, one needs bilayer graphene 18,19 , of which the energy dispersion becomes parabolic. Therefore, it is desirable to have Dirac materi-als with the properties of Dirac cones tunable by external experimental conditions. In this paper, we propose a new family of Dirac materials, including (RO) 2 (M X 2 ) 2 [20][21][22][23] and (AeF) 2 (M X 2 ) 2 24-26 (R = rare earth, Ae = Sr or Ba, M = Sb or Bi, X = S, Se or Te, O = oxygen and F = fluorine), with electrically tunable multiple Dirac cones.

II. ELECTRICALLY TUNABLE DIRAC CONES
To illustrate crystal structures of this family of materials, we may take (LaO) 2 (SbSe 2 ) 2 as an example, which possesses the tetragonal ZrCuSiAs-type structure 27 with space group P 4/nmm. As is shown in Fig. 1(a), (LaO) 2 (SbSe 2 ) 2 has a triple-layer (TL) structure with one (LaO) 2 layer sandwiched by two SbSe 2 layer in one unit cell. The (LaO) 2 layer is similar to that in LaOFeAs 28 , formed by a square lattice of O atoms, coordinated tetrahedrally with four neighboring La atoms. However, the SbSe 2 layer has a different structure from the FeAs layer. For one SbSe 2 layer, the Sb 1 (Sb ′ 1 ) atom and Se 2 (Se ′ 2 ) atom form a distorted checkerboard lattice while the Se 1 (Se ′ 1 ) atom lies between (LaO) 2 layer and the checkerboard SbSe layer. Each Sb atom can be viewed to have a distorted octahedral coordination with six Se atoms. Strong ionic bonding is formed between the central (LaO) 2 layer and two adjacent SbSe 2 layers within one unit cell while the chemical bonding between two adjacent SbSe 2 layers in neighboring unit cells is much weaker, dominated by van der Waals interaction 29 . The primitive lattice vectors are shown in Fig. 1(a), denoted as x , y , z . The lattice constant is denoted as a along the x and y directions and c along the z direction. The space group P 4/nmm is non-symmorphic and possess a glide operationĝ z = {m z | τ } which consists of a reflection in the xy-plane,m z : (x, y, z) → (x, y, −z), followed by a non-primitive translation τ = ( a 2 , a 2 , 0). The system has inversion symmetryÎ with the inversion center located at the center of two nonequivalent O atoms. In addition, there is a mirror symmetry along both the x and y directions, denoted asm x andm y . Based on the density-functional theory, we calculate the electronic band structure of (LaO) 2 (SbSe 2 ) 2 , which is found to be a narrow gap semiconductor with the bandgap near the X and Y points in the BZ, as shown in Fig. 1(b). Both conduction and valence bands are four-fold degenerate at X and split into two doubly degenerate states along the Γ-X line. An anti-crossing occurs between the conduction and valence bands along the Γ-X line, resulting in a bandgap around 20 meV. In addition, the band structure is found to be non-dispersive along the line Γ-Z and R-X, indicating that adjacent TLs are weakly coupled and the system is essentially two dimensional. Therefore, we will study a film configuration with one TL of (LaO) 2 (SbSe 2 ) 2 . It has recently been shown that a TL film is stable for one material (LaO) 2 (BiS 2 ) 2 in this family 29 . The configuration of a TL film is shown schematically in Fig. 2(a), with front and back gate voltages, which provide an asymmetric potential eU on the thin film. Fig. 2 The Fermi velocity is shown as a function of the gate voltage U . vA1 (vA2) indicates the velocity of Dirac cone A along the direction of kx (ky). from these two momenta. We emphasize that all the Dirac cones here are spin-resolved, different from that in graphene but similar to the case of TIs. Strikingly, the properties of Dirac cones, such as their positions, velocities and anisotropy, are tunable with gate voltages. For example, the velocity of Dirac fermions around X is shown as a function of the gate voltage U in Fig. 2(f), which can be tuned in the range (4.6 ∼ 6.7) × 10 5 m/s along the k x direction and (0.2 ∼ 4.0) × 10 5 m/s perpendicular to k x with the experimentally feasible gate voltages (the corresponding electric fields are in the range 0 ∼ 37 mV/Å 30,31 ). We notice that the velocity of linear dispersion has a small anisotropy for U < 200 mV, while the velocity anisotropy increases rapidly when U ≥ 200 mV, as clearly shown in the insets of Fig. 2(d) and (e). At the gate voltage U = 400 mV, the velocity along the Γ-X direction is around 7 × 10 5 m/s, about 4 times smaller than that in graphene, while that perpendicular to the Γ-X direction is one order smaller (2 × 10 4 m/s). To understand the existence and tunability of Dirac cones, we will next develop a low energy effective model for this system.

III. LOW-ENERGY EFFECTIVE MODEL
To develop a theoretical model, we first check the orbital natures of conduction and valence bands. The band structure with different atomic projections is shown in Fig. 3(a). The orbitals of La and O atoms lie far away from the Fermi energy due to the strong electron negativity and affinity. The outmost shells for both Sb (5s 2 5p 3 ) and Se (4s 2 4p 4 ) are p orbitals. The Se 1 and Se ′ 1 atoms are close to the (LaO) 2 layer and form strong bonds, which push their energy levels away from the Fermi energy. Thus, the bands near the Fermi energy are dominated by the p orbitals of Sb 1 (Sb ′ 1 ) and Se 2 (Se ′ 2 ) atoms in the SbSe checkerboard layers. Therefore, we focus on the bilayer SbSe checkerboard lattice and develop an atomic tight-binding (TB) model on this lattice with three p orbitals on each site. The details of lattice structures and the form of TB model are described in Supplementary Section B. By carefully choosing parameters, we can qualitatively reproduce energy dispersions near X (or Y) and four Dirac cones emerge along the Γ-X line after turning on an asymmetric potential between two layers, as shown in Fig. 3(b) and Supplementary Fig. S4.
Below, we will focus on the effective theory around X (or Y). We find that the conduction band minimum and valence band maximum are dominated by p x (p y ) orbitals of Se and Sb atoms around the X (Y) point. This orbital nature is confirmed by the ab initio calculation, as shown in Fig. S1 of Supplementary Section A. As described above, both conduction and valence bands have four-fold degeneracy at X, which originates from two spin states and two checkerboard SbSe layers. Therefore, we can construct a four band effective model around X on the basis |p x , σ, ξ µ for the conduction and valence bands respectively, where σ =↑ z , ↓ z denotes spin, ξ = c, v denotes conduction band and valence band, and µ = ± denotes the upper and lower layer. We choose the spin axis along the z direction. It should be emphasized that there is a strong hybridization between Se and Sb atoms for both the conduction and valence bands near the bandgap. Thus, we use ξ = c, v instead of Sb,Se to denote the basis. Using the Löwdin perturbation method 32 , we find the effective Hamiltonian takes the form around X up to the second order in k for the conduction or valence band, where ǫ ξ (k x , k y ) = d 0ξ + d 1ξ k 2 x + d 2ξ k 2 y , the Pauli matricesτ denote layer index andσ denote spin index. d iξ and f iξ are material dependent parameters, and can be extracted from perturbation procedure. In the Hamiltonian (1), the first term is not important. The second term is a spin-independent term, describing the hybridization between two layers. The third and fourth terms depend on spin and originate from the third order perturbation combining the interlayer hopping and SOC (see Supplementary Section C). We notice that these two terms take the familiar form of the Rashba type of Hamiltonian while the additionalτ z dependence indicates that the spin splitting is opposite for two layers. This layer dependent Rashba term is the origin of the layer dependent spin texture described below. The last term describes the asymmetric potential eU induced by gate voltages.
We next look at the symmetry properties of Hamiltonian (1). The present system has time reversal (TR) symmetryT and the space group symmetry P 4/nmm. The wavevector group of P 4/nmm at X contains glide symmetryĝ z , mirror symmetrym y and inversion sym-metryÎ. The effective Hamiltonian (1) can be derived from symmetry principles based on the above symmetries, as shown in Supplementary Section D. We find that the four-fold degeneracy at X can be determined from three symmetry operationsT ,ĝ z andm y . On the basis |σ, µ of the effective Hamiltonian (we neglect p x and ξ for short), the representation of symmetry operators is given by g z =σ zτx , m y = iσ y and T = +(−)iσ yτz K for conduction (valence) bands, where K is complex conjugate. Since [m y , H ξ ] = 0 at X, we may choose eigen states with definite mirror parity ofm y . Four eigen states can be written as |ψ µ,↑y(↓y) = 1 √ 2 (| ↑ z , µ + (−)i| ↓ z , µ ), which satisfym y |ψ µ,↑y(↓y) = +i(−i)|ψ µ,↑y(↓y) . Here ↑ y and ↓ y correspond to up and down spin along the y direction. Now let us consider howT andĝ z act on these four states |ψ µ,σy . BothT andĝ z change spin up ↑ y to spin down ↓ y . However, the obtained states operated bŷ T orĝ z are different since the layer index is preserved byT but changed byĝ z . In addition, the combination ofT andĝ z leads to the state with the same spin but different layer indices. Thus, these four eigen states can be related to each other byT ,ĝ z andTĝ z , so they must be degenerate at X. In general, the four-fold degeneracy at the X (Y) point is a direct consequence of the combination of TR symmetryT and the anti-commutation relation betweenĝ z andm y (m x ).
After understanding the degeneracy at X, we next consider the states away from X along the Γ-X line (k x = 0, k y = 0 in Hamiltonian (1)). In Fig. 2(b), one can see that four-fold degenerate states are split into two doubly degenerate states. According to Hamiltonian (1), both the hybridization term (f 1ξ term) and the SOC term (f 2ξ and f 3ξ terms) can contribute to this splitting. The remaining double degeneracy at a finite k x comes from the combined symmetryTÎ (Kramers' doublet due to both the TR symmetry and inversion symmetry). Alternatively, we can also understand it from the combination ofĝ z andm y . We can still choose two degenerate eigen-states to have definite mirror parity, denoted as m y |ψ ↑y(↓y) (k x ) = +i(−i)|ψ ↑y(↓y) (k x ) . It should be emphasized that the states |ψ ↑y(↓y) (k x ) should be a linear combination of the basis in different layers. |ψ ↑y (k x ) and |ψ ↓y (k x ) are still related to each other byĝ z which will also change the layer index. This means that if |ψ ↑y (k x ) mainly stays at the upper layer, |ψ ↓y (k x ) must be at the lower layer. Thus, spin and layer indices are related to each other, leading to the layer dependence of spin textures in this system. To show it more explicitly, we calculate spin polarization at a given momentum k for different layers based on the effective model (1).
As an example, we may consider spin texture of two degenerate eigen-states |ψ ξ α,+ (α = 1, 2) with eigen-energy E + = ǫ ξ + (f 2 1ξ + f 2 2ξ )k 2 x + f 2 3ξ k 2 y , which is given by where the layer dependent spin operator is defined as S i,µ = 2σ iτ z +µ 2 with i = x, y, z. From Eq. (2), we indeed find opposite spin textures for different layers once the parameter f 2ξ or f 3ξ , originating from SOC, is non-zero. Now let us turn on the asymmetric potential eU between two layers. Since spin states are locked to the layer indices, two spin states with opposite layer indices will be split accordingly. This picture has been utilized to explain the giant Rashba spin splitting in (LaO) 2 (BiS 2 ) 2 in the early study 29,33 . In our case, we need to consider spin splitting for both conduction and valence bands. It is found that the parameters f 2c and f 2v have opposite signs (f 2c > 0 and f 2v < 0) while f 3c and f 3v have the same sign (f 3c , f 3v < 0), leading to different spin textures for the conduction and valence bands. From the Supplementary Section E, for the case with U > 0, k x > 0 and k y = 0, the lowest conduction band carries spin up at the lower layer while the highest valence band is dominated by the states with spin down at the upper layer. Since these two states have opposite spin, as well as opposite mirror parities ofm y , there is no coupling between them. Thus, the Dirac cones due to the crossing points between the conduction and valence bands along the Γ-X line are protected by mirror symmetry. From the above analysis, we can see that the layer dependent spin texture is the underlying physical reason for the existence and tunability of Dirac cones. To further confirm our physical picture, we perform the first-principles calculations of layer dependent spin texture in this system, as shown in Fig. 3(c). The upper (lower) panel is the spin-resolved band projecting on the upper (lower) SbSe checkerboard layer at the gate voltage U = 100 mV. The opposite spin textures for the upper and lower SbSe layer are shown in the inset of Fig. 3(c), respectively.

IV. ELECTRICALLY TUNABLE QUANTUM ANOMALOUS HALL EFFECT
In the following, we will discuss the physical consequence of multiple Dirac cones in this system. We will consider the exchange coupling with magnetic moments and predict the quantum anomalous Hall (QAH) effect with a high Chern number in this system, which can be controlled by an external gate voltage. Magnetic moments can be introduced into this system by magnetic doping, which have been successfully used for TI materials, such as Sb 2 Te 3 family of materials, to realize the QAH effect 34 . Alternatively, one can also grow the TL film of (LaO) 2 (SbSe 2 ) 2 on top of ferromagnetic materials to induce exchange coupling by ferromagnetic proximity 35 . To theoretically study how exchange coupling affects our system, we construct a more realistic TB model using the maximum localized Wannier function method 36,37 , and then introduce exchange coupling phenomenologically. This method has been widely adopted to study the QAH effect in magnetically doped TIs 38,39 .
We start from the case with an asymmetric potential eU = 100 meV and no magnetization, in which eight Dirac cones exist in the whole BZ. The results obtained from the Wannier function method is consistent with the direct calculation from the first-principles methods. The exchange coupling is introduced into the calculation by the phenomenological Kondo-type of Hamiltonian H ex = J ex s · M , where J ex is the exchange coupling constant, s denotes electron spin and M is the average magnetization of the system. The Dirac cones are gapped by turning on magnetization (see Fig. S6 in Supplementary Section G). It is well-known that each gapped 2D Dirac cone (massive Dirac Hamiltonian) contributes half quantized Hall conductance ( e 2 2h or − e 2 2h ). Thus, the total Hall conductance is determined by the sign of the Hall conductance contribution from different Dirac cones. It turns out that all massive Dirac Hamiltonians take the same sign, leading to the total Hall conductance 4e 2 h or − 4e 2 h for the whole system. To see this more explicitly, we directly calculate edge states of the whole system in a ribbon configuration and plot the local density of states at one edge along the x direction. As shown in Fig. 4(a) and (b), there are in total four chiral edge states propagating along the same direction, with two appearing near theX point and the other two near theΓ point, whereΓ andX are the projection of Γ and X into the 1D edge. Therefore, the edge state picture is consistent with the analysis of bulk Dirac cones, revealing that the QAH state with the Hall conductance ± 4e 2 h can be realized in this system. Fig. 4(c) and (d) show bulk bandgap as a function of magnetization (J ex,Sb M z,Sb and J ex,Se M z,Se ) and asymmetric potential eU . Since a topological phase cannot be changed when the bulk bandgap remains open, the phase diagram can be determined by tracking the gap closing lines 40 . Two gapless lines are found in Fig. 4(c), indicating two topological phase transitions. For each metallic line, four equivalent Dirac cones (Dirac cones A or B in the inset of Fig. 2(c)) reverse their bandgap, leading to the change of Hall conductance by 4e 2 h . Therefore, we can determine the Hall conductance in each regime, as dictated by the regimes I, II and III in Fig. 4(c) with Hall conductance 4e 2 h , 0 and − 4e 2 h , respectively. We notice that the trivial insulating regime II with zero Hall conductance is quite small for (LaO) 2 (SbSe 2 ) 2 , but much larger for (SrF) 2 (SbSe 2 ) 2 (see Supplementary Section G). We emphasize that along the line a−b in Fig. 4(c), a topological phase transition from the Hall conductance 4e 2 h to 0 to − 4e 2 h can occur by tuning only gate voltages and fixing magnetization. This again reflects the electrical tunability of Dirac physics in this system.

V. DISCUSSION AND CONCLUSION
The physics discussed above for (LaO) 2 (SbSe 2 ) 2 can also be applied to other materials in this family. Since all low energy physics occurs in the SbSe 2 layer, the (LaO) 2 layer can also be replaced by other (RO) 2 layer where R is a rare earth atom. As shown in Supplementary Section A, the SbSe 2 layer can also be replaced by other SbX 2 layers (X = Te, S). For BiX 2 (X = S,Te,Se) 41,42 , an indirect bandgap occurs, ranging from 70 to 300 meV, and is difficult to be closed by tuning an external gate voltage. In addition, the (LaO) 2 layer can also be replaced by (AeF) 2 where Ae = Sr, Ba [24][25][26] . Experimentally, bulk (LaO) 2 (SbSe 2 ) 2 and (AeF) 2 (SbSe 2 ) 2 have been fabricated by the hightemperature ceramic method 20,24 . Although we focus on one TL film, we expect similar physics can be applied to multiple TLs. The electrical tunability of Dirac physics in this system indicates its potential application in various fields. For example, we have shown that a topological phase transition between the QAH states with different quantized Hall conductance can be achieved by controlling gate voltages, which will be useful for the experimental study of critical phenomena of topological phase transitions 34,43 . Moreover, superconductivity has been realized in (LaO) 2 (BiS 2 ) 2 [21][22][23]26 or in LaOFeAs, a well-studied unconventional superconductor with a similar crystal structure 44 . Therefore, it is possible to fabricate heterostructures combining (LaO) 2 (SbSe 2 ) 2 and these superconducting materials, which provide a new flatform to study the coexistence of Dirac physics and superconductivity.

VI. COMPUTATIONAL METHODS
All the first-principles calculations are based on the density-functional theory as implemented in the Vienna ab initio simulation package (VASP) 45,46 . The projector augmented wave method 47 is used, with a kinetic energy cutoff of 400 eV for the plane wave basis set. The generalized gradient approximation of the Perdew-Burke-Ernzerhof (PBE) type functional 48 is adopted to describe the exchange-correlation interaction. In the crystal structure, we take the experimental lattice constants: a = 4.13 A, and c = 14.17Å 20 . For the TL thin film, we build a slab model with a vacuum region of 14Å to decouple the consecutive slabs in the supercell approach. 13×13×5 and 13×13×1 Γ-centered k-point meshes are used in the bulk and slab calculations, respectively. The SOC is employed in all electronic structure calculations. We also perform electronic band calculations using the WIEN2K package 49 and reproduce energy dispersion obtained by VASP. The 18 × 18 tight binding Hamiltonian, which nicely captures the low energy physics around the Fermi level, is constructed using the maximal localized Wannier function method 36,37 . We choose p x , p y , p z states of Sb and Se atoms as the projection centers.

ACKNOWLEDGMENTS
We would like to thank Binghai Yan for useful discussions and Gang Yang for the help of the WIEN2K package. XYD and BFZ acknowledge the support from National Natural Science Foundation of China (Grant No.  11374173)