Enhanced valley splitting of WSe2 in twisted van der Waals WSe2/CrI3 heterostructures

Van der Waals (vdW) heterostructures composed of different two-dimensional (2D) materials offer an easily accessible way to combine properties of individual materials for applications. Owing to the discovery of a set of unanticipated physical phenomena, the twisted 2D vdW heterostructures have gained considerable attention recently. Here, we report enhanced valley splitting in twisted 2D vdW WSe2/CrI3 heterostructures. In particular, the splitting can be 1200% (or 5.18 meV) of the value for a non-twisted heterostructure. According to the k·p model, this value is equivalent to a ~20 T external magnetic field applied perpendicular to the 2D sheet. The thermodynamic stability of 2D vdW WSe2/CrI3 heterostructures, on the other hand, depends linearly on the interlayer twisting angle.


INTRODUCTION
Since the discovery of graphene 1 , some other atomically thin 2D materials, such as transition metal dichalcogenides (TMDs) 2-6 and 2D magnetic chromium triiodide (CrI 3 ) 7,8 , have gained attention because of their fascinating properties for applications. As a matter of fact, TMDs have a long historical standing as semiconductors with layered structures 9 . Early in 2010, Mak et al. 5 and Splendiani et al. 3 have revealed independently, that a strong photoluminescence emerges when MoS 2 crystal is thinned down to a monolayer, indicating an indirect to direct bandgap transition. Henceforth, there is a strong resurgence of interest in 2D TMDs, owing to the challenges and opportunities for applications in electronics 10 , optoelectronics 11 , and spintronics 2,12,13 .
Beyond individual materials, stacking different 2D materials into vdW heterostructures can give us more room to achieve improved functions 14,15 . Except for AA or AB stacking, 2D vdW twisted heterostructures 16 , in which one layer is rotated with respect to the other by an angle θ, are also interesting owing to their unexpectedly physical properties. For instance, Cao et al. 17 reported unconventional superconductivity emerging in 2D twisted bilayer graphene, and Wu et al. 18 showed topological insulators in 2D twisted TMDs homostructures. Similarly, Liao et al. 19 experimentally confirmed a~5 times enhancement of vertical conductivity in twisted MoS 2 /graphene heterostructures. These emerging twistronics merging spintronics could pave the way for future 2D material design and applications 20 .
In monolayer TMDs, a pair of degenerate but inequivalent energy valleys (K, K′) in momentum space, protected by time-reversal symmetry, become the third degree of freedom for information storage, other than charge and electron spin 21 . To break the timereversal symmetry, one of the effective approaches is applying an external perpendicular magnetic field (B) on TMDs, which can induce valley polarization through the Zeeman effect. Li et al. 22 obtained the magnitude of valley splitting of 0.12 meV T −1 in monolayer MoSe 2 , with B ranging from −10 to 10 T. Similarly, in monolayer WSe 2 , the valley splitting up to 0.25 meV T −1 is obtained 23,24 . However, a large external magnetic field is impractical for devices. Therefore, incorporating intrinsic 2D magnetic materials into vdW heterostructures, utilizing the large magnetic proximity effect (MPE), is an attractive alternative. It has been shown that MPE originating from 2D substrate can enhance the valley splitting 25,26 .
The recently fabricated magnetic monolayer CrI 3 7,27 with an out-of-plane easy axis is a suitable candidate to generate MPE when forming 2D vdW heterostructures with TMDs 28 . Zhong et al. 29 experimentally demonstrated that up to 3.5 meV valley splitting (which is equivalent to~10 T external B in single-layer WSe 2 ) has been achieved. Furthermore, this MPE can be tuned optically, which in turn alters the valley splitting of WSe 2 30 . Theoretical investigations also demonstrated that the valley degeneracy of WSe 2 can be lifted in a vdW WSe 2 /CrI 3 heterostructure because the time-reversal symmetry is broken by magnetic Cr ions [31][32][33] . Zhang et al. 31 showed a modulated energy splitting from 0.31 to 1.04 meV, which is determined by interlayer distance and/or atom arrangement. Hu et al. 32 confirmed that the valley splitting is sensitive to interlayer distance, increasing from 2.0 to 4.5 meV when the distance was decreased by 0.3 Å from its equilibrium value. In other words, in TMD/CrI 3 heterostructures, MPE is sensitive to interlayer vdW interaction. So far, however, the role θ played in 2D vdW WSe 2 /CrI 3 heterostructures is still unclear. How does the valley splitting depends on θ needs to be clarified.
In this work, we investigate a set of 2D vdW WSe 2 /CrI 3 heterostructures with different θ and stacking patterns with an emphasis on the magnetic properties, such as the magnetic moment and magnetic anisotropy energy (MAE). After that, we studied the electronic properties of 2D vdW WSe 2 /CrI 3 heterostructures by using a band unfolding method. The enhanced KK′ valley splitting of WSe 2 , having strong dependence on θ, is obtained. Finally, we analyzed the MPE effects through the difference of partial charge density, and deduced the equivalent B by a k·p model 34 .

Structural models and stability
Pristine WSe 2 has a (hexagonal) sandwiched structure, in which one W layer is sandwiched between two S layers. Pristine CrI 3 also has a sandwiched structure, with Cr layer being sandwiched between two I layers. We obtain the vdW WSe 2 /CrI 3 heterostructure by placing a WSe 2 layer on a CrI 3 sheet, with θ ranging between 0 and 30˚. To find stable configurations, we slide the WSe 2 layer relative to the CrI 3 layer. Figure 1a, b shows an example of the atomic structures with a 0˚twist angle. Here, we use ten points within a unit cell along with the zigzag (ZZ) and armchair (AC) directions, respectively. Figure 1c shows the energies along with the ZZ (upper) and AC (lower) directions, both suggesting that the starting structure is most stable. Here, the Cr atom is at the hexagonal center (HC) of WSe 2 , named Cr-HC hereinafter. In the following, we consider mainly the Cr-HC configurations for twisted WSe 2 /CrI 3 . However, we also consider other configurations, in particular, the less-stable Se-HC configurations where the Se atom is at the HC of Cr layer. The WSe 2 (CrI 3 ) monolayer belongs to the space group of P-6m2 (P-31m), which is reduced to P3 after stacking into the twisted heterostructures for either Cr-HC or Se-HC. The detailed atomic structures for the twisted Cr-HC structures can be found in the Supplementary Information.
For twisting, there are, in principle, many possibilities. Owing to the lattice constant difference (3.321 Å of WSe 2 vs 7.002 Å of CrI 3 ), most of them will result in huge supercells, making DFT calculations impossible. The lattice constant mismatch (δ) is defined as  3 36 are sensitive to the lattice constant. For this reason, in this work we mainly have performed two limited sets of calculations: first using the lattice parameter of WSe 2 and second using that of CrI 3 for the heterostructure. Furthermore, we investigated the strain effects on the thermodynamic stability and electronic properties for the heterostructures with the relaxed lattice parameter.   The stability of WSe 2 /CrI 3 is determined by their formation energy defined as where S is the cross-section of the supercell, E CrI3 , E WSe2 , and E CrI3/ WSe2 are total energies of respective systems. Table 1 below and Supplementary Fig. 1 show that twisted WSe 2 /CrI 3 have larger E f than that of non-twisted ones, indicating that the twisted heterostructures are more stable. Noteworthy, by using the same lattice parameters for individual and heterostructure systems, the strain energy owing to lattice mismatch have been excluded here. The energy differences between the Cr-HC and Se-HC configurations with the same twist angle θ range between 0.01 and 4.68 meV Å −2 , showing a high lubricity of 2D WSe 2 /CrI 3 similar to that of layered graphene 37 . Also noted is that the spin-orbit coupling (SOC) part of the formation energy also depends on θ.

Magnetic properties and valley splitting
For WSe 2 /CrI 3 with WSe 2 lattice constant, MAE decreases with increasing δ (from −5.14% of compressive strain to 3.36% of tensile strain), which is similar to monolayer CrI 3 38 , but it is independent of sliding between WSe 2 and CrI 3 . For WSe 2 /CrI 3 with a CrI 3 lattice constant, on the other hand, only a slight fluctuation of MAE is observed.
Next, we consider the valley splitting, schematically shown in Fig. 3. With SOC, both the (high-lying) valence band (VB) and (lowlying) conduction band (CB) of WSe 2 split into two subbands: one spin up and one spin down. In particular, the VB (CB) in the K valley will split into VB1 and VB2 (CB1 and CB2) and the VB′ (CB′) in the K′ valley will split into VB1′ and VB2′ (CB1′ and CB2′), denoted by dashed lines in Fig. 3b. When WSe 2 is put on a magnetic CrI 3 substrate, the spin-up bands (red) are shifted downward and the spin-down bands (blue) are shifted upward. The VB (CB) at K and K′ with different spins will split by Δ VB (Δ CB ), shown in Fig. 3b.
The effective KK′ valley splitting (Δ KK′ ) is defined as Note that to calculate the KK′ splitting, requires the calculation of energy levels at K (K′) in the unfolded Brillouin zone. This is done as follows: the wavefunction of the supercell (Ψ) can be mapped on to those of the primitive cell (Ψ k ) 39-41 by where N is the number of primitive cells contained in the supercell, k is the reduced wave vector inside the first Brillouin zone, is the integral multiple of lattice vectors a j ,T R is the translational operator for a translation R, and For a given region in WSe 2 /CrI 3 such as WSe 2 layer in the heterostructure, the relative weight (ρ) of layer projection in the range between Z 1 and Z 2 is given by as proposed by Chen and Weinert 42 . Figure 4 shows, for WSe 2 /CrI 3 with WSe 2 lattice parameter, the band structures for θ = 0°and 16.1°, (a, c) before and (b, d) after band unfolding. We see that although there is interaction between WSe 2 and CrI 3 , the unfolded band structures of WSe 2 look similar with a band gap of 1.24 ± 0.01 eV, close to that of monolayer WSe 2 . Table 2 lists the components for valley splitting in WSe 2 . For monolayer WSe 2 , both Δ VB (Δ CB ) and Δ KK′ are zero due to time reversal symmetry. When placed on CrI 3 , this symmetry is broken so valley splitting is generally observed. Consistent with previous calculations 31, 32 , at θ = 0°, |Δ KK′ | = 0.43 and 0.27 meV for Cr-HC and Se-HC, respectively. We find the enhanced valley splitting of twisted heterostructures, for either Cr-HC or Se-HC. Compared with θ = 0°, the magnitude can also be magnitude larger (e.g., 0.43 meV vs 5.18 meV). With the parameter relaxation, as presented in Supplementary Table 1, the valley splitting of twisted WSe 2 /CrI 3 heterostructure with Cr-HC is still larger than that of the non-twisted ones. We notice that the valley splitting without the twisting arises mainly from CB states, while for twisted ones, VB states also play an important role.
The enhanced valley splitting in heterostructures may be understood in terms of an MPE. This can be seen in Table 2, which shows that for both K and K′, E VB decreases, but E CB is nearly a constant, with an increasing θ. One may note that both the VB and CB states of WSe 2 are made of predominantly W d orbitals: d x 2 Ày 2 =xy for VB, and d z 2 for CB. Figure 5a-d shows the differences in the partial charge density Δρ of the valence band maximum (VBM) states at K and K′. An obvious θ-dependence of Δρ is observed near Cr atoms. Figure 5e shows the planar averaged Δρ along z, ΔρðzÞ. It reveals how the states of WSe 2 is affected by the proximity of Cr atoms inside CrI 3 . This analysis suggests that the VBM states are responsible for the enhanced valley splitting. Indeed, results in Table 2 support such a conclusion as it shows that valley splitting is mainly a result of a θ-dependent Δ VB . Figure 5f plots the sum of Δρ of Cr layer (ΣΔρ) as a function of Δ KK′ . A nearly linear dependence is observed. Since Δρ is a measure of proximity effect, Fig. 5f thus suggests that MPE is the reason for enhanced Δ KK′ . Figure 5f also shows that with an increase in valley splitting, the total magnetic moment of Cr also increases from 3.011 to 3.204 μ B . Hence, the enhanced magnetic    field is also a direct consequence of MPE via an enhanced valley splitting.
In the discussion above, we have chosen a lattice parameter, e.g., that of WSe 2 , to construct the supercells. It is desirable to consider the effect of strain. For this reason, Table 1 also shows the results when the lattice parameter is that of CrI 3 . In line with above results, Δ KK′ for the Cr-HC configuration shows an increase with θ, e.g., 1.97 meV (θ = 16.1°) < 3.69 meV (θ = 28.1°).
There are two important consequences due to strain: (1) a direct-to-indirect gap crossover between θ = 0 and 23.4°. To this end, recall that monolayer TMDs are direct gap semiconductors, in which both the VBM and conduction band minimum (CBM) are at the K (K′) valley of the Brillouin zone (see Fig. 3a) 3,5 . Previous experiments 43,44 and theory 35 also showed that the various physical properties of TMDs, including band gap and band edge positions in the Brillouin zone, can be tuned by applying an inplane strain. In our calculations, the strain in WSe 2 changes from 5.4% (tensile) at θ = 0°to −3.25% (compressive) at θ = 23.4°. In the former case, the CBM position is unchanged, but the VBM position moves from K (K′) to Γ. In the latter case, the VBM position is unchanged, but the CBM position moves to the Q (Q′) valley, which is about halfway between K (K′) and Γ, as shown in Fig. 3a.
(2) The valley splitting for the Se-HC configuration vanishes at θ = 16.1°. This is shown in Table 1 where for θ = 16.1°, MAE = 0.27 meV Cr −1 is positive with an in-plane easy axis. The corresponding in-plane magnetic field can neither generate a Zeeman splitting nor lift the valley degeneracy. However, in this case, it is possible to obtain an easy axis component perpendicular to the 2D sheet by applying an external electric field (E-field). In the absence of the E-field, a small dipole of −0.016 eÅ is obtained, which is in line with a 0.0037 eCr −1 transfer from WSe 2 to CrI 3 by a Bader analysis 45 . Figure 6 shows the results for both dipole moment and KK′ valley splitting. A good linear relationship of the dipole moment with applied E-field is observed, and Δ KK′ can be tuned obviously by applying the external electric field. Taking Efield = −0.3 V Å −1 as an example, we find that MAE changes from 0.27 to −0.56 meV Cr −1 , and Δ KK′ changed from 0 to −5.84 meV. This magnetoelectric effect here is reminiscent of what has been observed in bilayer CrI 3 46 , explaining the twist-induced enhancement of valley splitting in WSe 2 /CrI 3 heterostructures.
The k·p model To estimate the magnitude of the equivalent B and understand the Zeeman effect induced by CrI 3 , we use a k·p model 34 , in which the interaction energy is divided into spin and orbital contributions: where g 0 = 2 is the Landé factor, μ B ¼ e h=2m 0 is the Bohr magneton, m 0 is the electron mass, S = σ=2 is the spin operator with σ being the Pauli matrices, and L is the orbital angular momentum operator. As mentioned earlier, in the K (K′) valley the CBM state mainly consists of W d z 2 orbital with L Z = 0, whereas the VBM state mainly consists of W d x 2 Ày 2 =xy orbital with L Z = 2. From the DFT computation, the spin projection around CB and VB band edges remains mostly out of plane with little in-plane tilting as shown in Supplementary Table 2. Therefore, the Rashba interaction is neglected in the present k·p model. The value P i;j VBM i jH B jCBM j , where i and j run over the four levels near the Fermi level as shown in Fig. 3b, can be directly calculated via DFT, from which we deduce the equivalent B. For WSe 2 /CrI 3 heterostructures with the Cr-HC configuration and WSe 2 lattice parameter, we obtain B = 2.1 T for θ = 0°, 7.4 T for θ = 16.1°, and 20.9 T for θ = 23.4°. It shows that not only the magnetic field strength is sufficiently strong but also twisting can be an effective way to amplify the MPE at interfaces.

DISCUSSION
In summary, using first-principles calculation, we perform a systematic study on the structural stability, electronic and magnetic properties of 2D vdW WSe 2 /CrI 3 heterostructures. We show that, compared to the non-twisted structure, the twisted ones are more stable, as evidenced by their higher interfacial binding energies. Our study also reveals an MAE of the heterostructures, which is controlled by strains of the CrI 3 layer. The MPE can be understood in terms of the charge density difference of the VBM states at K and K′. More importantly, in twisted heterostructures, there is an order of magnitude enhancement of K-K′ valley splitting, which can also be tuned by applying an external electric field. With the help of a k·p model, we determine the equivalent B field due to twisting as the origin of an enhanced MPE.

METHODS
First-principles calculations were performed using the Vienna ab initio simulation package (VASP) based on the density functional theory (DFT) 47 . We employed the Perdew-Burke-Ernzerhof (PBE) 48 functional, within the projector augmented wave (PAW) 49 approach, for exchange-correlation potential and energy. The plane-wave energy cutoff was set at 500 eV. vdW interactions between CrI 3 and WSe 2 were included by employing Grimme's semiempirical DFT-D3 scheme 50 . To avoid interaction between periodic images in the supercell approach, we used a vacuum space larger than 30 Å. The atomic structures were fully relaxed until the Feynman-Hellman force on each atom is less than 0.02 eV Å −1 . The electronic self-consistent convergence criterion was set at 10 −7 eV. The K-point mesh 51 in the Brillouin zone was sampled with the density of 0.03 Å. To determine the magnetic anisotropy, the spin-orbit coupling (SOC) was included in our DFT calculations. In addition, we performed electronic structure calculations with GGA + U methods 52 described by Dudarev, in which the on-site Coulomb parameter U and the moderate exchange parameter J were set to 2.7 and 0.7 eV 53 , respectively. The k-projection method 42 for interfaces modeled by supercells within the framework of the first-principles method was used to obtain the unfolded electronic band structures.

DATA AVAILABILITY
The data that support this work are available in the article and Supplementary Information file. Further raw data are available from the corresponding author (J.Z.) on reasonable request.