Electric field induced giant valley polarization in two dimensional ferromagnetic WSe2/CrSnSe3 heterostructure

Valleytronics is receiving extensive research efforts. Thus, we investigated the electric field-induced valley polarization in the WSe2/CrSnSe3 heterostructures by varying the stacking order. The heterostructure shows indirect band gaps of 270 and 330 meV in the two most stable structures. The WSe2/CrSnSe3 heterostructure displays a ferromagnetic ground state with out-of-plane anisotropy (0.02 meV) in one stable stacking (S-1) while a small in-plane anisotropy (−0.01 meV) is found in other stacking (S-2). The Curie temperature is slightly enhanced to 73 K compared to the monolayer CrSnSe3. We have found the valley splitting of 4 meV in S-1 whereas it became 9 meV in the S-2 system. The valley splitting is further enhanced if an electric field is applied from CrSnSe3 to the WSe2 layer whereas it is suppressed in the reversed electric field. Particularly, the S-2 structure shows a giant valley splitting of 67 meV at an electric field of 0.6 V Å−1. We attribute this electric field-dependency to the dipolar effect. Overall, we propose that the WSe2/CrSnSe3 heterostructure can be a potential structure for obtaining a giant valley splitting.


INTRODUCTION
Study on two-dimensional (2D) materials has become one of the most fascinating issues in condensed matter physics and materials sciences because many intriguing physical quantities have been found in 2D materials. In these 2D materials, along with the traditional charge and spin degree of freedom physical quantities, another physical quantity namely valley state has been introduced as another degree of freedom which can be utilized for future quantum information storage devices 1 . The valleytronics (valleydependent electronics) is receiving tremendous research interest because its quantum information storage may exceed the existing charge or spin-based technology. Regarding the valleytronics, numerous exciting phenomena such as valley-dependent circular dichroism 2 , valley-polarized electroluminescence 3 , valley Hall effect 4 , spin-valley locking effect 5 , valley excitons and intervalley quantum coherence are reported 6 . Among many two-dimensional (2D) materials for valleytronics, mostly the transition metal dichalcogenide (TMD) systems have been explored. For instance, the TMDs such as MoS 2 , MoSe 2 , WS 2 , and WSe 2 have broken inversion symmetry in the odd number of layer structures whereas the inversion symmetry is preserved in the even number of layer structure [7][8][9][10] . Usually, the odd number of layer TMD materials have a strong spin-orbit coupling (SOC) effect in the valence bands due to the lack of inversion symmetry, and this SOC effect is mostly observed for Group 6 of TMD materials 11 . Therefore, these materials are treated as potential candidates for valley related studies. Due to the broken inversion symmetry, the energy levels are splitted at K + and Kvalley but still both the valleys are energetically degenerate since the timereversal symmetry is still preserved 12 . This degeneracy will be lifted in the presence of a magnetic field due to the broken time reversal-symmetry. This feature is so-called the valley polarization, and the valley polarization is crucial in the valley related devices because we can selectively produce the carriers from inequivalent valleys by applying an electric field. So far, several approaches are proposed to lift the valley polarization by applying external magnetic field 13,14 , doping magnetic ions 15 , creating vacancies 16 , or incorporating the magnetic proximity effect in a heterostructure [17][18][19] . However, the valley splitting under an external magnetic field is not practical for device applications because a huge magnetic field is required to split the valley state. For instance, only a small valley splitting of 0.1-0.2 meV T −1 was generated in previous reports 13,20,21 . Alternatively, magnetic impurity doping has been proposed to investigate the valley physics, but the magnetic impurity doping may induce metallic band structure 15 . Besides, the breaking of valley degeneracy was also reported by using an intense circularly polarized light through the Stark effect 22,23 , but the restricted lifetime of carriers limit their practical applications for information storage 2 . Among the several approaches, the van der Waals heterostructure can be a promising method because of its stability and easy controllability 24 . For instance, Zhang et al. reported the stacking effect on the valley polarization in the WSe 2 /CrI 3 heterostructure 25 . Similarly, the effect of twisting the CrI 3 monolayer and gate voltage in MoSe 2 /CrI 3 and WSe 2 /CrI 3 heterostructure was also studied by Zollner et al. 17 . Similarly, we also studied the valley splitting and electric field tuning of valley splitting in graphene/CrI 3 heterostructures 26 . So far, the reported intrinsic valley polarization (valley splitting) in CrI 3 /TMD van der Waals heterostructure is rather weak ranging from about 1~3 meV. Thus, it is highly desirable to fabricate a structure that displays a giant valley splitting. To this end, we will investigate the valley splitting in the WSe 2 /CrSnSe 3 heterostructure, and propose the electric field induced giant valley splitting.

CrSnSe 3 and WSe 2 monolayers
We first explored the magnetic properties of monolayer CrSnSe 3 before we investigate the valley polarization due to the magnetic proximity effect. Thus, the 2D CrSnSe 3 monolayer structure is constructed from the experimentally reported 2D CrGeTe 3 (CGT) crystal structure by replacing the Ge and Te atoms with Sn and Se. After replacing these atoms, we performed the structure relaxation without any constraints. Figure 1(a) shows the top view of the optimized crystal structure of CrSnSe 3 monolayer. Here, the purple, pink, and gray sphere represents the Cr, Sn, and Se atoms. Note that the lattice constant of 2D CGT is 6.920 Å. However, the optimized lattice constant of CrSnSe 3 monolayer was 6.608 Å, and this suppression was mainly due to the difference in the atomic radii of Ge, Sn, Se, and Te. The optimized CrSnSe 3 monolayer belongs to the space group of P-31M (#162) so that every Cr atom is bonded with six Se atoms. The distance between two Cr atoms was 3.815 Å and the interatomic distance between Cr and Se atoms was 2.643 Å forming a distorted octahedron. On the other hand, each Sn atom was bonded with three Se atoms at a distance of 2.60 Å in the upper layer and one Sn atom in the lower layer at 2.704 Å. To check the stability of the monolayer CrSnSe 3 structure, we calculated the formation energy by using the following relation where E total , µ Cr , µ Sn , and µ Se are the total energy of the CrSnSe 3 monolayer, the chemical potentials of one Cr, Sn, and Se atom. To calculate the chemical potentials, we considered the energy of a single element from the body-centered cubic crystal structure for bulk Cr, face-centered cubic for Sn, and hexagonal crystal for Se. We found that E total = −48.10 eV for CrSnSe 3 monolayer while µ Cr , µ Sn , and µ Se were −14.05 eV, −7.94 eV, and −20.39 eV. Using the above relation, we obtained the formation energy of −0.57 eV per atom. We also examined the dynamical stability by calculating the phonon dispersion spectrum. For phonon calculations, the crystal structure was relaxed until the forces on every atom are less than 0.005 eV Å −1 . The phonon band structure along the high symmetry directions is shown in Supplementary Fig. 1a. Since no imaginary frequencies appear, we confirm that the 2D CrSnSe 3 monolayer is dynamically stable. Along with the dynamical stability, we also investigated the thermal stability of monolayer CrSnSe 3 using molecular dynamic (MD) simulations. Here, we considered two ensembles in the MD simulation: (i) NVE ensemble with fixed number of atoms (N), volume (V) and energy (E) and (ii) NVT ensemble with fixed number of atoms (N), volume (V) and temperature (T). In Supplementary Fig. 1b, c, we present the top and side views of the snapshot of the obtained atomic structure at 8 ps at 400 K using both NVE and NVT ensembles. The structural integrity remains intact even after 8 ps at 400 K. Also, we show the evolution of temperature (K) and the corresponding energy (eV) at 400 K using both NVE and NVT ensembles in Supplementary Fig.  2a-d. Here, both NVE and NVT ensembles showed almost the same results. This feature may indicate the thermal stability of the 2D CrSnSe 3 monolayer at least up to 400 K. To investigate the magnetic ground state, we calculated the total energy difference between ferromagnetic (FM) and the antiferromagnetic (AFM) spin configurations. The FM ground state was obtained with an energy difference (E ex = E AFM − E FM ) of 91 meV per cell, and the total magnetic moment was 6 µ B per unit cell. The interatomic distance between two Cr atoms was 3.815 Å and the bonding angle of Cr-Se-Cr became 92.4˚. This may indicate that the FM state between Cr atoms is mediated by Se atom through superexchange interaction suggested by Goodenough-Kanamori-Anderson (GKA) rule [27][28][29] . The local magnetic moment of Cr, Sn, and Se atoms in the pristine CrSnSe 3 layer were 3.38, 0.04, and −0.15 µ B . To understand the electronic band structure nature of the CrSnSe 3 ferromagnetic layer, we present the electronic band structure of pristine CrSnSe 3 monolayer in Fig.  1b. Here, the blue lines represent the majority spin (spin-up) bands while the red lines are for the minority spin (spin-down) bands. The 2D monolayer CrSnSe 3 had an indirect bandgap of 0.72 eV, and the band gap was slightly decreased to 0.63 eV after including the SOC. The conduction band minimum (CBM) was originated from the majority spin band, and the Cr d-orbitals with a small hybridization with p-orbitals of Sn and Se atoms mostly formed the CBM band. Meantime, the valence band maximum (VBM) was built from the minority spin band of p-orbital of the Sn and Se atoms. We also calculated the magnetocrystalline anisotropy. Here, we performed a non-collinear total energy calculation including the spin-orbit coupling along [010] directions. The pristine CrSnSe 3 monolayer had a uniaxial magnetic anisotropy (perpendicular to the 2D sheet, i.e., along [001] direction) of 0.03 meV per unit cell. According to the SOC matrix analysis, both Cr and Sn had uniaxial anisotropies of 0.06 and 0.01 meV per atom while four Se atoms had small in-plane anisotropies of −0.03 meV per atom. Note that two other Se atoms had no contribution to the magnetic anisotropy. Consequently, we obtained a small perpendicular magnetic anisotropy of 0.03 meV per cell in the pristine CrSnSe 3 layer. Besides, we calculated the Curie temperature (T C ) using Metropolis Monte Carlo (MC) simulations. In this simulation, the temperaturedependent magnetization curve was calculated using 10000 equilibration steps and 10000 averaging steps by using the VAMPIRE software package. We considered a supercell of 50 × 50 with periodic boundary conditions in the MC simulations to reduce the finite-size effects. The obtained value of the Curie temperature was 61 K, and this is slightly larger than found in the monolayer CrI 3 (45 K) 30 . On the other hand, we fully optimized the WSe 2 monolayer without intralayer van der Waal (vdW) interaction and obtained lattice constants of a = b = 3.32 Å. These values agree with the experimental report 31 . Figure 1c shows the two dimensional honeycomb lattice of WSe 2 , and Fig. 1d shows the projected band structure after including the SOC. Here, we consider the perpendicular magnetization to the monolayer surface. A direct band gap of 1.25 eV is obtained at K + and K − valleys, and this agrees with the previously reported value 32 . The VBM consists of the d x 2 -y 2 and d xy orbitals while the d z 2 orbital of the W atom form the CBM band. Since the WSe 2 monolayer has the time-reversal symmetry, we obtained the same SOC induced splitting at K + and K − valleys in both valence and conduction bands, and the opposite spin state is found at both valleys due to the broken inversion symmetry. This implies that the spin and valley indices are locked at the band edges. On the other hand, Hsu et al. reported an indirect band gap nature in the WSe 2 monolayer structure 33 . Indeed, the WSe 2 monolayer has three sublayers. So, we further investigated the band structure including intralayer vdW interaction to elucidate the conflicting results. With including the vdW, we obtained lattice parameters of a = b = 3.29 Å and an indirect band gap of 1.28 eV. The obtained band structure is presented in Supplementary Fig. 3. Overall, we found that the lattice parameters were dependent on the existence of vdW interaction, and the direct or indirect band gap nature of the monolayer WSe 2 monolayer seems sensitive to the change in the lattice parameters. WSe 2 /CrSnSe 3 heterostructure We now focus on the WSe 2 /CrSnSe 3 heterostructure. As described above, the lattice constant of the CrSnSe 3 layer was 6.608 Å while it was 3.32 Å for WSe 2 monolayer. Thus, we employed the (2 × 2) supercell of WSe 2 monolayer. Since the lattice mismatch is only 0.5%, we use the lattice constant of 6.608 Å in the heterostructure. Note that the vdW interaction in the WSe 2 /CrSnSe 3 heterostructure is included for structure optimization. Thus, the disagreement in the band gap nature of the WSe 2 monolayer mentioned above does not affect the calculated physical properties of the WSe 2 / CrSnSe 3 heterostructure. To find the most stable stacking order, we considered six different stackings denoted by S-1, S-2, S-3, S-4, S-5, and S-6 in Fig. 2. In S-1, one W was located between two Cr atoms whereas one W and one Se were on top of two Cr atoms in S-2. In S-3, one Se in the WSe 2 layer was placed between two Cr atoms while one W was on top of the upper Se of CrSnSe 3 layer in S-4. In S-5, one W was on top of the lower Se of CrSnSe 3 layer and one Se of WSe 2 was on top of Se of CrSnSe 3 layer in S-6 stacking. Through the total energy calculations, we found the most stable state in S-1 stacking. Nonetheless, the total energy difference between S-1 and S-2 was only 8 meV per unit cell. Meantime, the total energy difference with S-3 stacking was 106 meV, and also the energy differences with other structures (S-4, S-5, and S-6) were in the range of 20-24 meV. In the real experimental situations, the total energy difference less than 10 meV may be overcome. Thus, we focus on the S-1 and S-2 structures in our study. The interlayer distance (d) in S-1 stacking was 3.23 Å while it is 3.25 Å in S-2 stacking. We also calculated the binding energy (E b ) of the heterostructures using the following relation where E tot , E CrSnSe3 , and E WSe2 are the total energies of the optimized heterostructure, isolated monolayer CrSnSe 3 , and isolated WSe 2 layer while N is the total number of atoms in the heterostructure. We found that E tot = −142. calculated binding energies were −45.07 meV per atom in S-1 and −44.69 meV per atom in S-2 structures. These energies are comparable to the binding energies reported in CrI 3 /WSe 2 heterostructure 25 . Like in the monolayer CrSnSe 3 structure, we found a total magnetic moment of 6 µ B per cell in both the S-1 and S-2 stackings, and this is due to the weak van der Waal interaction between WSe 2 and CrSnSe 3 layers. We obtained a perpendicular magnetic anisotropy (PMA) of 0.02 meV per cell in S-1 stacking whereas the S-2 stacking displayed an in-plane magnetic anisotropy of −0.01 meV per cell. Since the in-plane magnetic anisotropy energy is quite weak, it will be possible to switch from in-plane to perpendicular direction by applying a small external magnetic field. In the heterostructure, we found that the total energy difference between FM and AFM was slightly enhanced from 91 meV in the monolayer to 108 meV in the bilayer heterostructure. This results in the enhancement of the Curie temperature to 73 K.
Since the central issue in this work is to investigate the valley polarization. Thus, we calculate the band structures after including the SOC. Note that the vdW interaction was included in the band structure calculation. The direction of magnetization is assumed to be perpendicular to the surface. Figure 3a, b shows the projected band structures of S-1 and S-2 stacked systems. The blue and red colors represent the up and down spin components of the band. In S-1 stacking, we obtained an indirect band gap of 270 meV while the S-2 displayed an indirect band gap of 330 meV. Besides, we have presented the projected band structure of S-1 and S-2 structures in Supplementary Fig. 4. We obtained that the VBM was mainly originated from the WSe 2 layer whereas the CBM was originated from the CrSnSe 3 layer. This indicates that the WSe 2 / CrSnSe 3 heterostructure belongs to type-II heterostructure 34 . The band gap suppression in the heterostructure will be discussed later. Figure 3c, d shows the unfolded band structures in S-1 and S-2 stacked systems projected on WSe 2 monolayers. The valley feature of the WSe 2 monolayer was still preserved in the heterostructure. From this band structure, we calculated the valley splitting by using the relation.
where E c and E v is the energy of the conduction and valence bands at K + and K − valleys. The black arrow represents the direction of the magnetization of the 2D magnetic layer. In comparison with the bare WSe 2 layer, the degeneracy at K + and K − valleys was broken, and we also found a stacking order dependency. For instance, the S-1 stacking had a valley splitting of around 4 meV while the S-2 system showed a valley splitting of 9 meV. Note that the valley splitting in both systems (S1 and S2) of the WSe 2 /CrSnSe 3 heterostructure is larger than that found in the previous reports of the WSe 2 /CrI 3 heterostructure (~3.5 meV) 35 . To reveal the origin of stacking dependent valley splitting, we calculated the average charge density in the x-y plane along the zdirection and the charge density difference. The calculated results are presented in Fig. 4a, b and insets. We found no meaningful charge transfer although the charge redistribution was observed in the heterostructure. Thus, the electrostatic potential gradient appears as displayed in Fig. 4c, d, and this resulted in a vertical electric field along the z-direction. Therefore, we found electric dipole moments of 0.20 and 0.26 Debye in S-1 and S-2 stackings pointing from the CrSnSe 3 to the WSe 2 layer. In the heterostructure, the CBM (originated from the CrSnSe 3 layer) and the VBM (originated from the WSe 2 layer) were dragged towards the Fermi level due to this internal dipolar field, and the band gaps in heterostructures were suppressed compared with their individual component. Since both inversion and time-reversal symmetries are broken in the heterostructure due to the magnetic layer, the valley splitting in the WSe 2 may occur due to the proximity effect created by the magnetic CrSnSe 3 layer. Note that the valley degeneracy at K + and K − has recovered if the magnetization direction lies in the in-plane direction. Thus, the valley splitting in the WSe 2 layer can only be observed if the effective magnetic field from the CrSnSe 3 layer is perpendicular to the WSe 2 layer. Besides, we found that the magnitude of valley splitting had a stacking order dependency. This stacking dependency in valley splitting  can be attributed to the different nearest neighboring (NN) distances between W and Cr in two stackings. For instance, the NN distance between W and Cr in S-1 stacking was 6.83 Å however, in S-2 this distance was 6.58 Å. A small distance between W and Cr may result in a strong coupling and a relatively large valley splitting in S-2 stacking compared with the S-1 stacked system 25 .
To enhance the valley splitting, we applied an external electric field. As discussed earlier, the dipole moment was pointing from the CrSnSe 3 to the WSe 2 layer at zero electric field in both S-1 and S-2 structures. Under the external electric field, we can expect charge redistribution so that the dipole moment will have an electric field dependency. Figure 5a shows the dipole moment under the external electric field. The negative electric field means that the electric field is applied from the CrSnSe 3 to the WSe 2 layer while the positive electric field is the reversed direction. The dipole moment was further enhanced in the negative electric field while the direction of the dipole moment was changed under the positive electric field. The enhancement in dipole moment with increasing the electric field strength is due to the more inclination of electronic charge toward the CrSnSe 3 or WSe 2 layers. As an illustration, we presented the charge distribution affected by the electric field in S-1 stacking, and the inset of Fig. 5a shows the result. The electronic charge was more inclined towards the CrSnSe 3 layer under the negative electric field while it was more inclined toward the WSe 2 under the positive electric field. Overall, the band gap was suppressed in the negative electric field while it was enhanced in the positive electric field as shown in Table 1. Figure 5b, c shows the projected band structure including SOC in [001] direction at −0.6 V Å −1 field in S-1 and S-2 structures. Note that the projected band structures at other electric field strengths in S-1 and S-2 stackings are presented in Supplementary Fig. 5 and Supplementary Fig. 6. Compared with the band structure at zero electric field displayed in Fig. 3a, b, both CBM and VBM moved toward the Fermi level due to the enhancement of the dipolar field, and the band gap was decreased to 3 and 52 meV in S-1 and S-2 structures. Figure 5d shows the magnified band structure near the K + edge for both S-1 and S-2 stackings. This change in band structures can also affect the valley splitting. Figure 5e shows electric field induced valley splitting in both S-1 and S-2 stackings. The magnitude of valley splitting was strongly dependent on the electric field strength and direction. Particularly, we obtained a giant valley splitting of 67 meV at a negative electric field of −0.6 V Å −1 while the valley splitting was decreased to 5 meV in the positive electric field of +0.6 V Å −1 in S-2 structure. We attribute this to the dipolar effect. To illustrate this, we calculated the charge density difference in S-2 structure at the electric field of −0.4 V Å −1 and +0.4 V Å −1 . Figure 6a-c shows the top view of the charge density difference at zero, −0.4 V Å −1 , and +0.4 V Å −1 electric fields. As shown in Fig. 6a at zero field, the charge distribution had a spatial dependency along with the K + and K − directions of the Brillouin zone. Thus we expect a different local dipolar field at K + and K − valleys, and this gives rise to the valley splitting in pristine systems. This anisotropic behavior in the charge density difference was further enhanced under the negative electric field as shown in Fig. 6b. Thus, we obtained increased valley splitting. In contrast, the anisotropic feature was suppressed under the positive electric field in the K + and K − direction as shown in Fig. 6c, and this resulted in the reduced valley splitting. We also calculated the magnetocrystalline anisotropy under the external electric field. In the S-1 system, the magnetic anisotropy was insensitive to the electric field. In S-2 stacking the magnetic anisotropy remained insensitive under the positive electric field. However, we observed the spin reorientation from out-of-plane anisotropy to small in-plane anisotropy under a negative electric field as shown in Table 1.

DISCUSSION
In summary, we explored the structural stability, electronic and magnetic ground state of CrSnSe 3    the S-1 structure displayed a valley splitting of 4 meV whereas it became 9 meV in the S-2 system. The valley splitting was further enhanced if the electric field was applied from the CrSnSe 3 to the WSe 2 layer, and it was suppressed in the reversed electric field. Particularly, the S-2 structure showed a giant valley splitting of 67 meV at an electric field of −0.6 V Å −1 . We attribute this electric field dependency to the dipolar effect. Overall, our findings may imply that the WSe 2 /CrSnSe 3 heterostructure can be a potential structure for obtaining a giant valley splitting.

First-principles calculations details
We performed the spin-polarized density functional calculations using the Vienna ab initio simulation package (VASP) 36,37 . The generalized gradient approximation (GGA) of Perdew-Burke-Ernzerhof (PBE) was used as an exchange and correlation functional in both CrSnSe 3 monolayer and the WSe 2 /CrSnSe 3 heterostructures 38 . A Hubbard type onsite Coulomb interaction was employed for the correlation effect among d-electrons of Cr atoms. The Dudarev et al. 39 formalism with a value of U eff = 3 eV was selected for d-orbitals of Cr atom. Besides, Shen et al. reported the spin polarized band structure of the monolayer CrSnSe 3 using the U = 34 eV 40 . We included the van der Waals correction DFT-D3 41 in the heterostructure. The vacuum thickness of 20 Å was employed to avoid the artificial interaction between the periodic images along the z-axis. A plane wave energy cutoff of 600 eV was used for all the calculations. The lattice constants and internal coordinates of the 2D CrSnSe 3 monolayer were fully optimized with energy and force convergence of less than 10 −5 eV and 0.01 V Å −1 by using the conjugate gradient method. The first Brillouin zone was sampled by using a Γ centered k-point mesh of 11 × 11 × 1 for structure relaxation and electronic properties. The optimized lattice constant of the CrSnSe 3 monolayer was applied for the WSe 2 /CrSnSe 3 heterostructures because the lattice mismatch was quite weak (less than 0.5%). However, the internal coordinates were fully relaxed in all systems.

Phonon dispersion and ab initio molecular dynamics
The dynamical stability of the monolayer CrSnSe 3 was explored by calculating the phonon dispersion using a finite difference approach within PHONOPY 42 . For phonon dispersion curve, we employed the 3 × 3 × 1 supercell with a k-point mesh of 5 × 5 × 1, and the force criterion for the ionic step was set to 10 −8 V Å −1 . We also investigated the thermal stability by using ab initio molecular dynamics (MD) 43 where the Nose thermostat was used in both NVE and NVT ensembles 44 at 400 K temperature lasted for 8 ps with a time step of 0.5 fs. A supercell of 2 × 2 × 1 and a k-point mesh of 3 × 3 × 1 is used in MD simulation.

Magnetic anisotropy and Curie temperature
The total energy method was used to calculate the magnetocrystalline anisotropy. The magnetocrystalline anisotropy was considered as a difference in the total energy when the spins are aligned along [001], [010] and [100] directions in presence of spin-orbit-coupling 45 . A denser kpoint mesh of 21 × 21 × 1 was used in all the studied systems for calculating the magnetic anisotropy and the convergence of magnetic anisotropy energy was carefully checked. We also investigated the effect of the vertical electric field (ranging from 0.2 to 0.6) on the magnetocrystalline anisotropy and valley splitting in different stackings. The Curie temperatures were calculated for both CrSnSe 3 monolayer and bilayer heterostructure using Metropolis Monte Carlo (MC) simulations 46,47 .

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.