Magnetic properties of bilayer VI3: Role of trigonal crystal field and electric-field tuning
Thi Phuong Thao Nguyen, Kunihiko Yamauchi, Tamio Oguchi, Danila Amoroso, Silvia Picozzi
MMagnetic properties of bilayer VI :Role of trigonal crystal field and electric-field tuning Thi Phuong Thao Nguyen, Kunihiko Yamauchi, ∗ and Tamio Oguchi Institute of Scientific and Industrial Research, Osaka University,8-1 Mihogaoka, Ibaraki, Osaka, 567-0047, Japan
Danila Amoroso and Silvia Picozzi
Consiglio Nazionale delle Ricerche (CNR-SPIN),Unit`a di Ricerca presso Terzi c/o Universit`a “G. D’Annunzio”, 66100 Chiety, Italy (Dated: February 9, 2021)The magnetic properties of two-dimensional VI bilayer are the focus of our first-principles analy-sis, highlighting the role of trigonal crystal-field effects and carried out in comparison with the CrI prototypical case, where the effects are absent. In VI bilayers, the empty a g state - consistent withthe observed trigonal distortion - is found to play a crucial role in both stabilizing the insulating stateand in determining the inter-layer magnetic interaction. Indeed, an analysis based on maximally-localized Wannier functions allows to evaluate the interlayer exchange interactions in two differentVI stackings (labelled AB and AB’), to interpret the results in terms of virtual-hopping mechanism,and to highlight the strongest hopping channels underlying the magnetic interlayer coupling. Uponapplication of electric fields perpendicular to the slab, we find that the magnetic ground-state inthe AB’ stacking can be switched from antiferromagnetic to ferromagnetic, suggesting VI bilayeras an appealing candidate for electric-field-driven miniaturized spintronic devices. I. INTRODUCTION
Boosted by the experimental discovery of intrinsicmagnetism in atomically thin layers of CrI , Cr Ge Te and Fe GeTe , two-dimensional (2D) van der Waals(vdW) magnets have recently received an increasing at-tention. Interestingly, the control of 2D magnetism infew atomic layers is enabled by external electric field or by electron-hole doping, making them particularlyappealing for potential spintronic applications.Among those materials, VI , belonging to the family oftransition-metal (M) trihalides MX (X = Cl, B, and I)with honeycomb arrangement of the metal cations (simi-lar to the most studied CrI ), has recently emergedas a potential 2D ferromagnet . It is known sincemore than 30 years ago that, in the bulk form, VI be-comes ferromagnetic below a Curie temperature of T c (cid:39)
55 K, similar to CrI (with T c (cid:39)
68 K) . Conversely,the structural properties are still under debate. Exper-imental characterizations of the crystal structure havein fact reported that VI undergoes a structural phasetransition around 79 K, changing its symmetry acrossstill unclear phases: the high-temperature (HT) crystalstructure was proposed to be either trigonal P c orrhombohedral R , or monoclinic C m structure ; thelow-temperature (LT) crystal structure was proposed tobe either C c or R ¯3 structure . Optical and electricaltransport measurements have clearly showed bulk VI tobe a semiconductor with an optical band gap of ∼ . . However, from the theoretical point of view, theunderstanding and modeling of the electronic propertiesare the focus of present debate. Some studies have infact reported that bulk VI is a Mott-insulator with abandgap of about 1 eV , whereas others have claimeda half-metallic character . In a thin-film limit, to the best of our knowledge, noexperimental studies have been reported for atomic lay-ers of VI . On the theoretical side, current studies arecontroversial with respect to the electronic properties,in analogy to the situation for the corresponding bulkphase. In particular, by analysing electronic propertiesin the VI monolayer, in Ref. a Mott-insulator groundstate is proposed, reported to be lower in energy than thehalf-metallic state (by ∼ have proposed orbital-ordered phases accompanying the lattice distortion, whilethe authors of Ref. ascribed the gap opening to com-bined effects of spin-orbit coupling (SOC) and Hubbard U correlations. On the other hand, a consensus is reachedwith respect to the magnetic properties of a VI singlelayer: current theoretical characterizations report a fer-romagnetic (FM) exchange coupling between first V-siteneighbors, i.e. FM intra-layer coupling. The inter-layermagnetic stability has also been investigated for bilayerVI . In particular in Ref. it is reported that the inter-layer magnetic stability is sensitive to the layer stack-ing, in line with previous works on bilayer CrI . Forthe sake of completeness, we mention that in Ref. it isclaimed bilayer VI to show a stacking-independent ferro-magnetic ground state, but considering the half-metallicstate rather than the Mott-insulating one.A deep understanding of the VI electronic structureis needed to interpret the related magnetic properties,both in the monolayer and in the bilayer case. In thisstudy, we therefore focus on the crucial role of trigonalcrystal field effects in determining the VI insulating be-haviour and the related magnetic properties. As a coun-terexample, we consider the prototypical 2D magnet, i.e. CrI (where trigonal crystal field effects are absent), and a r X i v : . [ c ond - m a t . m t r l - s c i ] F e b we carry out a one-to-one comparison between VI andCrI . In particular, we first focus on the monolayer anddiscuss crystal field effects (Section III A) and density ofstates (Section III B). Then we move to the the magneticproperties of bilayer halides, by considering two differentstacking arrangements. We concentrate on the interlayerexchange coupling, interpreting it in terms of virtual hop-ping mechanisms, highlighting the most efficient hoppingpaths (Section IV A) and addressing the effects of an ex-ternal electric field in tuning the magnetic stability (Sec-tion IV B). Finally, we draw our conclusions in SectionV. II. COMPUTATIONAL METHOD
Density-functional theory (DFT) calculations wereperformed using the VASP code within the general-ized gradient approximation (GGA) . The van derWaals (vdW) interactions were included for bilayer struc-ture calculations. The rotationally invariant GGA + U method was employed to account for correlation effects .On-site Coulomb interaction for transition-metal 3 d or-bitals was considered with an effective U of 2.0 eV .Semiconducting and metallic states are initialized via thedensity matrix within the GGA+ U scheme . Brillouinzone integrations were performed using a k -point grid of6 × × × × k -points mesh. Elec-tric fields are applied perpendicularly to the surface bysaw-tooth-like potential with dipole correction . Themaximally-localized Wannier functions (MLWFs) werecalculated by using WANNIER90 tool interfaced withthe VASP code. III. RESULTS FOR MONOLAYER VI A. Crystal-field effects
First we focus on the difference between VI and CrI monolayers, as far as crystal field effects are concerned.In particular, both VI and CrI show the magnetic atomcoordinated to six I atoms, forming edge-sharing octahe-dra and resulting in octahedral crystalline electric field(CEF) splitting of the d orbitals into the two-fold e g andthree-fold t g states. The t g state can be further splitinto a doublet e (cid:48) g and a singlet a g state to reduce theband energy, i.e. a Jahn-Teller (JT) effect in trigonalsymmetry . This effect leads to the trigonal latticedistortion, but such JT-distortion is not remarkable inthis system. The I-Cr-I bond angle in CrI is almost 90 ◦ ( i.e. a cubic octahedron), while the I-V-I angle in VI ap-proaches 89 ◦ . As such, it exhibits a trigonal distortion -elongation along the z axis - still preserving spatial inver-sion symmetry [see Fig. 1(a) and (b)]. Fig. 1(c) and (d) Cr (b)(d) O h e g t g z elongation D (a)(c) O h e’ g a g e g e g t g JT effect V e g e g a g e’ g e’ g (e) zx y (f) t g e g z’ y’ x’ CrVI
FIG. 1. Distortion of the crystal structure in (a) VI and(b) CrI . The distortion does not change the bond length butalters the bond angle, leading to a trigonal elongation alongthe z -direction. (c) and (d) show the crystal field splitting of d levels in V 3 d and Cr 3 d , respectively. Five 3 d Wannierfunctions reflect the trigonal CEF states in monolayer VI (e)and cubic CEF states in monolayer CrI (f). The isosurfacelevels of the Wannier functions were set at 1 . a − / (yellow)and − . a − / (blue), where a is the Bohr radius. show the difference in CEF splitting for VI and CrI re-spectively: as schematically represented, the JT-inducedsplitting allows for the band gap opening in VI by half-filling the majority e (cid:48) g channel and leaving the a g emptyin the case of V d . On the other hand, CrI is unaffectedboth because of the Cr d valence and the almost cubicCEF.Within the global Cartesian { xyz } coordinate system,Fig. 1(e), the a g and e (cid:48) g states are written in the form | a g (cid:105) = 3 z − r (cid:12)(cid:12) e (cid:48) g (cid:11) = 1 √ (cid:16) √ (cid:0) x − y (cid:1) − zx (cid:17)(cid:12)(cid:12) e (cid:48) g (cid:11) = 1 √ (cid:16) √ xy + yz (cid:17) (1)where the z axis is parallel to the out-of-layer direction( i.e. perpendicular to the slab). According to the dif-ferent local crystal field effects, we projected the Blochfunctions onto the local M I octahedral coordinate sys-tem { x (cid:48) y (cid:48) z (cid:48) } for CrI (with basis axes directed along theCr-I bonds), and onto the Cartesian { xyz } system forVI [see Fig. 1(e) and (f)]After the projection and maximally localization pro-cess, the Wannier functions converged into localized or-bitals, as shown in Fig. 1 (e) and (f): orbitals shapes arein agreement with the e g -( a g , e (cid:48) g ) states splitting inducedby the trigonal CEF in VI , and the e g - t g states inducedby the cubic CEF in CrI . In particular, for the latter,it is possible to recognize the (cid:12)(cid:12) z − r (cid:11) orbital shape forthe a g state (occupying the empty space at the centerof the I triangle and pointing along the z -direction) andmixed shapes from the ( (cid:12)(cid:12) x − y (cid:11) , | zx (cid:105) ) and ( | xy (cid:105) , | yz (cid:105) )orbitals for the e (cid:48) g and e (cid:48) g states respectively, accordingto Eq.(1). B. Density of states
In line with previous works, our DFT calculationson monolayer VI converged to two different elec-tronic states, corresponding to half-metallic and Mott-insulator states : the two V d electrons occupy t g states as a g e (cid:48) g and e (cid:48) g a g , respectively. We found thatthe insulating state (with a direct band gap of ∼ U = 2 . et al. proposed that the insulatingground state of monolayer VI is stabilized by spin-orbitcoupling (SOC) splitting rather than by CEF splitting.here we remark that the crystal field splitting is energeti-cally dominant over the SOC effect, as usually seen in 3 d transition-metal compounds . In fact, when includ-ing SOC in our band structure calculations, we observedthat it affects the width and energy bands related to theI- p states, but it does not significantly change the V 3 d band structure. Moreover, as we focus here on the e (cid:48) g a g state, eventual SOC induced splitting would not affectthe mechanisms and conclusions presented in this study.The partial density of states (DOS) for monolayer CrI and VI are shown in Fig. 2; DOS are resolved for eachMLWF state, clearly showing the D d and O h CEF split-ting for VI and CrI respectively, therefore validatingour basis functions choice for the Wannier projection asexplained in the following.We note that the Cr d partial DOS, projected on thetrigonal D d basic set, so as to allow a direct comparisonwith the V d case, shows an overlap in the energy rangeof the a g and e (cid:48) g ; this reflects the absence of splittinginduced by the trigonal CEF, in line with the O h localoctahedral symmetry of CrI . In particular, Cr- d orbitalstates are split into empty e g and occupied t g stateswith a gap of about 0.9 eV in the majority spin chan-nel (up-spin states); the minority spin channel (down-spin states), unoccupied for both orbital types, does notdisplay any relevant splitting. Such a different behav-ior between the majority and minority spin channels canbe ascribed to the pd hybridization: the up-spin d -statesstrongly hybridize with I- p states located at the top of thevalence band, causing the large CEF splitting supported DO S ( s t a t e s / e V / a t o m ) Energy (eV)-4 -2 (b) Cr d e g t (a) V d e g e’ g a FIG. 2. The partial DOS (within U = 2.0 eV) projected onto(a) V-3 d - D d and (b) Cr-3 d - O h CEF states, respectively, viaWannier functions. Blue (red) color represents for majorityspin (minority spin). Fermi level is set as energy origin. by the e g - p bonding-antibonding splitting; the down-spin d -states are higher in energy, i.e. away from I- p levels,thus not showing any significant CEF splitting.In VI , the spin-up channel of the V d states are clearlysplit into e g , a g , e (cid:48) g trigonal CEF states. In particular,the e (cid:48) g is the lowest energy state with a broad distribu-tion due to the pd hybridization below the Fermi level,similar to CrI . On the other hand, the CEF splittingshows a different behaviour of the empty d states: the a g state becomes the lowest energy state, while the e (cid:48) g statestill lies in the same energy region as the minority Cr- t g state. This is related to the fact that the e (cid:48) g state hasmore bonding character with surrounding I p state thanthe a g state (compare the orbital shapes in Fig. 1(e));therefore, the pd hybridization shifts unoccupied e (cid:48) g levelup and occupied I- p level down. IV. RESULTS FOR BILAYER VI A. Electronic properties and inter-layer magneticstability
Let us now consider the case of bilayer VI , showing asimilar atomic structure to its bulk counterpart. In par-ticular, we studied two structures corresponding to R ¯3and C c phases in bulk VI . The difference between b (a) (b)(c) (d)TopBottom b bcc aa b c a J ’ AB stacking AB’ stacking J J J ’ J ’ J ’ (e) (f) FIG. 3. (a, b) Top views and (c, d) side views of atomicstructure in bilayer VI in AB and AB (cid:48) stacking. The redand orange hexagons represents honeycomb structure of Vatoms in top and bottom layers, and gray balls represent Iatoms. The black arrow indicates the vector which connectsequivalent atoms located in two layers and shows how the toplayer is sliding with respect to the bottom layer. Inter-layerexchange coupling J ij in bilayer VI for (e) AB and (f) AB (cid:48) stacking. these two structures is related to the different stackingof the two VI single layers. As shown in Fig. 3 (a) and(c), the AB-stacking is characterized by the top layershowing one V atom sitting above the hexagon center ofthe bottom layer, similar to bilayer graphene. Indeed,when comparing equivalent atoms in two layers, the toplayer is horizontally shifted from the bottom layers by(2 a + b ) / a and b de-note the lattice vectors. This structure has R ¯3 symmetry.On the other hand, the AB (cid:48) stacking, with C m space-group symmetry, is characterized by the shift of the toplayer by − ( a + b ) /
2, as shown in Fig. 3 (b) and (d). Weemployed the experimental values for the in-plane lat-tice constants of both materials: in bulk VI , we used a = b = 6 .
84 ˚A and for CrI , we used a = b = 6 .
87 ˚A .A 20 ˚A-thick vacuum is contained in the supercell for 2Dslab simulation. TABLE I. Relative total energy (meV/f.u) for inter-layer FMand AFM spin configurations in bilayer VI and CrI in ABand AB (cid:48) stacking. The lowest energy state is highlighted.AB AB (cid:48) VI FM CrI FM AFM 14.7 36.6
First, we remark that the robustness of the intra-layerFM spin ordering is demonstrated by calculating the en-ergy difference between FM and N´eel type antiferromag-netic AFM spin configurations in VI monolayer, i.e. ∆ E = E AFM − E FM = 12.8 meV/f.u. The V magneticmoment was calculated as 2.16 µ B . We then address themagnetic properties of bilayers, by calculating the totalenergy between inter-layer FM and AFM orders for thetwo stackings in bilayer VI and compared the resultswith those obtained in CrI . In bilayer VI , the inter-layer FM and AFM spin configurations are very close inenergy; nevertheless, the FM order is favored in the AB-stacking, while the AFM order is favored in AB (cid:48) stack-ing, as reported in Table I. Differently, in bilayer CrI ,the FM order is favored in both AB and AB’ stackingpattern. Noteworthy, the inter-layer FM order in AB (cid:48) stacking is only slightly more stable than the AFM or-der, the energy differences being rather sensitive to usedon-site Coulomb U values, therefore not allowing a directcomparison with previous works on bilayer CrI . In anycase, according to the energy differences reported in Ta-ble I, the magnetic stability in bilayers VI and in CrI AB (cid:48) -stacking results to be weak. As such, it may leadto an easy control of the magnetism by either externalelectric fields or electrostatic doping.To understand the magnetic stability, we evaluated themagnetic exchange interactions between V atoms by fit-ting total energies calculated in AB and AB (cid:48) stacking tothe Heisenberg Hamiltonian. Here we assume the Heisen-berg Hamiltonian, H = (cid:88) (cid:104) i,j (cid:105) J ij s i · s j , (2)where J ij are the isotropic Heisenberg coupling constantbetween spin sites i and j and s i is the unit vector point-ing to the direction of the spin at site i . A parallel spin(FM) configuration is favored when J <
J > J (cid:107) ), we thus considered inter-layercouplings ( J and J in AB stacking; J (cid:48) , J (cid:48) , J (cid:48) and J (cid:48) inAB (cid:48) stacking) as schematically illustrated in Fig. 3 (e)and (f); associated atomic pairs distances are reportedin Table II. In particular, we performed calculations toestimate J (cid:48) , J (cid:48) and J (cid:48) in a 2 × × . This method allowsto consider one specific pair of spins and remove the back-ground interactions, therefore allowing the calculation ofthe inter-layer magnetic exchange coupling constants ofinterest.In Table II we report the estimated exchange couplingconstants for bilayer VI and CrI .For VI , the intra-layer exchange coupling favors par-allel spin state, while the inter-layer coupling eventuallyfavors parallel spin state in AB-stacking and anti-parallelspin states in AB’-stacking, (cfr Table I). In closer detail,in AB stacking, J favours anti-parallel coupling (0 . J favours parallel coupling ( − .
24 meV).
TABLE II. Number of equivalent bonds per unit cell N , bond distance between transion-metal sites d , and calculated exchangecoupling constants J ij in AB and AB (cid:48) stacking for bilayer VI and CrI .AB stacking AB (cid:48) stacking J (cid:107) J J J (cid:107) J (cid:48) J (cid:48) J (cid:48) J (cid:48) VI N d (˚A) 3.95 6.66 7.74 3.95 7.04 7.05 8.07 8.95 J ij (meV) -3.20 0.81 -0.24 -4.46 0.10 0.21 0.04 -0.04CrI N d (˚A) 3.95 6.57 7.68 3.95 7.00 7.02 8.03 8.92 J ij (meV) -7.03 -0.82 -0.69 -8.11 -0.18 -0.23 -0.29 0.25 Since there is one J bond and nine J bonds per unitcell, overall the ferromagnetic configuration is more sta-ble. In AB (cid:48) stacking, both J (cid:48) and J (cid:48) favour anti-parallelcoupling (0 . .
21 meV), thus dominantly con-tributing to the inter-layer AFM coupling stability. ForCrI , the intra-layer and inter-layer exchange couplingbasically favor parallel spin states in both AB stackingand AB (cid:48) stacking.In order to shed light on the microscopic mechanismbehind the stacking-dependent magnetic couplings, werecall the “virtual hopping” idea based on the Hubbardmodel . In particular, we discuss the results in termsof the virtual inter-layer hopping of e g - t g states sup-ported by the inter-layer M-I-I-M super-exchange effect.In the weak hopping limit, the inter-site hopping can betreated as a perturbation to the ground state in whichmagnetic ordering does not affect the energy. When thehopping process is allowed between occupied and unoc-cupied states, it in turn contributes to the ground stateenergy through the second-order contribution as the ef-fective exchange energy J eff = 2 t /U with hopping inte-gral t and Coulomb repulsion U , the process being called“virtual hopping”. If we consider the direct hopping be-tween occupied and unoccupied 3 d states at the transi-tion metal sites, the parallel-spin configuration is favoredif the hopping is strong between majority- and majority-spin states; on the other hand, the anti-parallel spin con-figuration is favored if the hopping between majority-and minority-spin states is strong. In order to discussthe virtual hopping process, we extracted the hoppingparameters by using a MLWF basis set, as illustratedin Fig. 1 (e) and (f). Note that the Wannier functionsare centered at V and Cr sites and spreading the tailto I sites, so that our virtual hopping process implicitlyincludes the pd hybridization process. The same con-cept can be found in Anderson’s original work on super-exchange interaction .Figure 4 shows the inter-layer hopping paths with thecorresponding MLWFs which are responsible for the ex-change energy in bilayer VI and CrI . Here we selectthree types of inter-layer exchange couplings: first neigh-bor and second neighbor ( J and J ) interactions in ABstacking; first neighbor ( J (cid:48) ) interaction in AB’ stacking.The calculated hopping integrals corresponding to theseexchange couplings are shown in Table. III.In AB-stacking bilayer VI , the trigonal CEF levels AB stacking AB’ stacking e ’ g2 e g1 e ’ g2 e ’ g2 VI a (a) 1 st neighbor (b) 2 nd neighbor (c) 1 st neighbor -1.6 -3.10.0 a (d) 1 st neighbor (e) 2 nd neighbor (f) 1 st neighbor d zx d z d xy d x -y d yz CrI d x -y FIG. 4. MLWFs relevant to inter-layer exchange couplingin bilayer (a-c) VI and (d-f) CrI . ↑ and ↓ denote the ma-jority and minority spin state, respectively. The arrows showthe electron hopping from an occupied orbital state to anunoccupied orbital state; the dashed and solid lines denotethe parallel-spin and anti-parallel-spin configurations, respec-tively. Values of the hopping integrals (meV) are also shownnearby the arrows. Isosurface levels were set at ± . a − / for (a-c) and ± . a − / for (d-f). and the two-electron occupation make J positive (anti-parallel-spin-favored). The hopping between e (cid:48)↑ g and a ↑ g states are calculated to be negligible ( t ∼ e (cid:48)↑ g and e ↓ g states issizable ( t =1.0 meV), which may be responsible for theanti-parallel-spin-favored exchange interaction. In AB-stacking bilayer CrI , in contrast, the negative (parallel-spin-favored) J can be explained by a sizeable hop-ping between occupied d ↑ zx state and unoccupied d ↑ x − y state ( t =0.6 meV). This is consistent with previousworks, claiming that the e g - t g hopping leads to theFM coupling . As shown in Fig. 4(d), the diagonallyelongated lobes of d zx and d x − y orbitals show a paththrough Cr-I-I-Cr sites with d ↑ zx - p – p - d ↑ x − y hybridiza-tion, where the first pd hybridization shows π -like andthe second shows σ -like bonding. The second neighborinteraction J is negative both for VI and CrI . This canbe explained by a large hopping integral between e (cid:48)↑ g and e ↑ g and that between d ↑ xy and d ↑ z states. In Fig. 4(b),we can recognize a e (cid:48)↑ g - p - e ↑ g hybridization, where the pd hybridization shows σ bonding. For the CrI case, a sim-ilar picture holds (cfr Fig. 4(e)), where the e (cid:48) g orbital isreplaced by d xy orbital. DO S ( s t a t e s / e V ) Energy (eV) Δ E ( m e V / V a t o m ) FM(a) (b)FMAB stacking AB’ stacking-0.10.00.11.01.11.21.30.0 0.1 0.2Electric field (V/Å)(c)(d) e ’ g a e ’ g a a tt a e g AFM0.0 0.1 0.2-2 20 1-1 0.2
FIG. 5. The energy difference between the FM and AFMordering for (a) AB and (b) AB (cid:48) stacking as a function ofan external electric field. The positive value of ∆E meansFM is favored and negative means AMF is favored. d -orbitalprojected density of state for top (solid filled) and bottom(solid line) layer of bilayer VI with AFM AB (cid:48) stacking inAFM ordering (c) without and (d) with external electric field.Black arrow presents for the energy shifted by applied electricfield. Vertical dash line denote the Fermi energy. The AB’-stacked bilayer VI shows an interesting in-terplay between the atomic arrangement and the specificcharacter of the a g orbital of V, related to the VI trig-onal CEF: since a vanadium atom in the bottom layer islocated right under an iodine atom in the upper layer, theV- a g orbital strongly overlaps with the I- p z orbital andform the σ bonding (Fig. 4(c)). This makes the e (cid:48) g - a g hopping relevant; in particular, the e (cid:48)↑ g - p – a ↓ g hopping isvery strong ( t = − . J (cid:48) positive. In AB’stacked bilayers, the inter-layer exchange interactions areweaker than those in AB stacked bilayers, since severalpossible hoppings between multiple orbital states tend to cancel each other due to the atomic arrangement. Inevident contrast, AB’-stacked bilayer CrI shows resultssimilar to the AB-stacked case: the d ↑ yz - p – p - d ↑ x − y hy-bridization path, where the former (latter) pd hybridiza-tion shows π ( σ ) bonding, makes J (cid:48) negative (Fig. 4(f)). B. Electric field control of magnetic stability
Finally, we discuss the effect of an applied electric fieldon the magnetic stability in bilayer VI . This is indeedrelevant, since a magnetic phase transition upon electric-field application has been reported in bilayer CrI .The energy difference between inter-layer AFM and FMstates in AB and AB’ stacked with applied electric fieldsis shown in Fig. 5. In both stacking cases, an appliedelectric field promotes the FM ordering. Remarkably, inthe AB (cid:48) stacking, the ground state switches from AFM toFM ordering when the electric fields exceed a thresholdvalue of ∼ d orbital stateof top and bottom layers in AB (cid:48) stacked bilayer VI isshown in Fig. 5 (c) and (d). Without electric field, theDOS relative to the top and bottom layer lie in the sameenergy range. As discussed above, there is a competitionbetween parallel-spin hopping and anti-parallel-spin hop-ping in determining the first-neighbor exchange coupling J (cid:48) . Since the energy difference between e (cid:48)↑ g and a ↓ g state(∆ ε ↑↓ = 2 . e (cid:48)↑ g and a ↑ g state (∆ ε ↑↑ =1.4 eV), one may thinkthat a parallel-spin configuration is favored. However, J (cid:48) is found to be slightly AFM-favored. This is because theanti-parallel spin hopping ( t ↑↓ = − . t ↑↑ = − . J eff ∝ t / ∆ ε . Upon electric fields, the a g orbital stateof the top layer is shifted up, while it is shifted down inthe bottom layer (cfr Fig. 5). The band gap becomes nar-rower due to the shift of DOS and in turn decreases thedifference of orbital energy levels, while the hopping in-tegral is not significantly affected. Overall, this increasesthe tendency toward FM stability, eventually switchingthe favored magnetic configuration from AFM to FM. V. CONCLUSIONS
By means of first-principles calculations, we investi-gated the magnetic stability in bilayer VI and comparedour results with the corresponding well-studied case ofCrI . In particular, the magnetic exchange interactionshave been analyzed by evaluating the hopping integralsbetween MLWFs projected onto 3 d orbital states at Vand Cr sites. We found out that the trigonal crystal field TABLE III. Hopping integrals calculated by MLWF basis set between occupied and unoccupied d orbital states in parallel-( t ↑↑ ) or anti-parallel ( t ↑↓ ) spin configurations. Three types of hopping integrals, t , t , and t (cid:48) , corresponding with inter-layerexchange couplings J , J , and J (cid:48) are listed. ∆ ε labels the difference between two eigenenergy for the MLWFs relevant to thehopping process, corresponding to ∆ ε ↑↑ and ∆ ε ↑↓ in the main text. The dominant hopping values relevant to the exchangecouplings and those illustrated in Fig. 4 are highlighted in bold.VI Hopping t ↑↑ Hopping t ↑↓ a g - a g e g - e g e (cid:48) g - a g e (cid:48) g - e g a g - a g e g - e g e (cid:48) g - a g e (cid:48) g - e g e (cid:48) g - e (cid:48) g ∆ ε (eV) 0 0 1.4 1.5 1.4 2.1 2.9 3.4 3.5 t (meV) -3.6 2.1 -0.7 -4.6 2.7 0.0 1.0 0.6 t (meV) -0.7 1.1 0.3 -1.6 -0.5 -0.9 -0.6 2.6 -1.0 t (cid:48) (meV) 3.6 1.0 -1.1 1.6 4.9 1.1 -3.1 Hopping t ↑↑ Hopping t ↑↓ e g - e g t g - e g t g - t g e g - e g t g - e g t g - t g ∆ ε (eV) 0 1.5 0 2.7 4.3 4.6 t (meV) 0.9 t (meV) 0.6 t (cid:48) (meV) 2.1 associated to the relevant JT distortion within a single-layer of VI (and absent in CrI ), plays an importantrole for the inter-layer magnetic exchange interaction.The t g orbital states are in fact split into e (cid:48) g and a g states; the latter shows the typical lobe shape pointingalong the out-of-plane direction and the strong hoppingbetween bottom-layer a g and top-layer e (cid:48) g states deter-mine the antiferromagnetic inter-layer coupling in AB’stacking bilayer VI . Nevertheless, since the hoppingsfavoring parallel-spin and anti-parallel-spin configurationare highly competing, the application of electric fields al-lows the switching of the inter-layer magnetic orderingfrom AFM to FM, paving the way to spintronic applica-tions of VI -based 2D magnets. ACKNOWLEDGMENTS
This work was supported by “Center for SpintronicsResearch Network”, Osaka University. The numerical computation was performed on the Supercomputing Fa-cilities at the Institute for Solid State Physics, Universityof Tokyo. KY also acknowledges Center for Computa-tional Materials Science, Institute for Materials Research,Tohoku University for the use of MASAMUNE-IMR(Project No.20K0045). SP and DA acknowledge financialsupport from the Italian Ministry for Research and Edu-cation through the Progetto Internazionale NanoscienceFoundry and Fine Analysis (NFFA-MIUR) facility andthrough the PRIN- 2017 project “TWEET: Towards Fer-roelectricity in two dimensions” (IT-MIUR Grant No.2017YCTB59).
REFERENCES ∗ current address: Elements Strategy Initiative for Catalystsand Batteries (ESICB), Kyoto University, Kyoto 615-8245,Japan B. Huang, G. Clark, E. Navarro-Moratalla, D. R. Klein,R. Cheng, K. L. Seyler, D. Zhong, E. Schmidgall, M. A.McGuire, D. H. Cobden, W. Yao, D. Xiao, P. Jarillo-Herrero, and X. Xu, Nature , 270 (2017). C. Gong, L. Li, Z. Li, H. Ji, A. Stern, Y. Xia, T. Cao,W. Bao, C. Wang, Y. Wang, Z. Q. Qiu, R. J. Cava, S. G.Louie, J. Xia, and X. Zhang, Nature , 265 (2017). Y. Deng, Y. Yu, Y. Song, J. Zhang, N. Z. Wang, Z. Sun,Y. Yi, Y. Z. Wu, S. Wu, J. Zhu, J. Wang, X. H. Chen, andY. Zhang, Nature , 94 (2018). B. Huang, G. Clark, D. R. Klein, D. MacNeill, E. Navarro-Moratalla, K. L. Seyler, N. Wilson, M. A. McGuire, D. H.Cobden, D. Xiao, W. Yao, P. Jarillo-Herrero, and X. Xu,Nature Nanotechnology , 544 (2018). S. Jiang, J. Shan, and K. F. Mak, Nature Materials ,406 (2018). E. S. Morell, A. Le´on, R. H. Miwa, and P. Vargas, 2DMaterials , 025020 (2019). S. Jiang, L. Li, Z. Wang, K. F. Mak, and J. Shan, NatureNanotechnology , 549 (2018). M. A. McGuire, H. Dixit, V. R. Cooper, and B. C. Sales,
Chemistry of Materials , Chemistry of Materials , 612(2015). N. Sivadas, S. Okamoto, X. Xu, C. J. Fennie, and D. Xiao,
Nano Letters , Nano Letters , 7658 (2018). D. Soriano, C. Cardoso, and J. Fern´andez-Rossier, SolidState Communications , 113662 (2019). P. Jiang, C. Wang, D. Chen, Z. Zhong, Z. Yuan, Z.-Y. Lu,and W. Ji, Phys. Rev. B , 144401 (2019). J. Kim, K.-W. Kim, B. Kim, C.-J. Kang, D. Shin, S.-H.Lee, B.-C. Min, and N. Park,
Nano Letters , Nano Letters , 929 (2020). D.-H. Kim, K. Kim, K.-T. Ko, J. H. Seo, J. S. Kim, T.-H.Jang, Y. Kim, J.-Y. Kim, S.-W. Cheong, and J.-H. Park,Phys. Rev. Lett. , 207201 (2019). S. Son, M. J. Coak, N. Lee, J. Kim, T. Y. Kim, H. Hami-dov, H. Cho, C. Liu, D. M. Jarvis, P. A. C. Brown, J. H.Kim, C.-H. Park, D. I. Khomskii, S. S. Saxena, and J.-G.Park, Phys. Rev. B , 041402(R) (2019). S. Tian, J.-F. Zhang, C. Li, T. Ying, S. Li, X. Zhang,K. Liu, and H. Lei,
Journal of the American ChemicalSociety , Journal of the American Chemical Society ,5326 (2019). Y.-P. Wang and M.-Q. Long, Phys. Rev. B , 024411(2020). C. Huang, F. Wu, S. Yu, P. Jena, and E. Kan, Phys.Chem. Chem. Phys. , 512 (2020). K. Yang, F. Fan, H. Wang, D. I. Khomskii, and H. Wu,Phys. Rev. B , 100402(R) (2020). J. He, S. Ma, P. Lyu, and P. Nachtigall, J. Mater. Chem.C , 2518 (2016). T. Kong, K. Stolze, E. I. Timmons, J. Tao, D. Ni, S. Guo,Z. Yang, R. Prozorov, and R. J. Cava, Advanced Materials , 1808074 (2019). J. A. Wilson, C. Maule, P. Strange, and J. N. Tothill,Journal of Physics C: Solid State Physics , 4159 (1987). C. Long, T. Wang, H. Jin, H. Wang, and Y. Dai,
The Jour-nal of Physical Chemistry Letters , The Journal of PhysicalChemistry Letters , 2158 (2020). G. Kresse and J. Furthm¨uller, Phys. Rev. B , 11169(1996). J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.Lett. , 3865 (1996). S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J.Humphreys, and A. P. Sutton, Phys. Rev. B , 1505(1998). A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, Phys.Rev. B , R5467 (1995). J. Neugebauer and M. Scheffler, Phys. Rev. B , 16067(1992). A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Vander-bilt, and N. Marzari, Computer Physics Communications , 685 (2008). W. B. Wu, D. J. Huang, J. Okamoto, A. Tanaka, H.-J.Lin, F. C. Chou, A. Fujimori, and C. T. Chen, Phys. Rev.Lett. , 146402 (2005). D. I. Khomskii,
Transition Metal Compounds (CambridgeUniversity Press, 2014). J. B. Goodenough, Phys. Rev. , 466 (1968). P. Bruno, Phys. Rev. B , 865 (1989). H. J. Xiang, E. J. Kan, S.-H. Wei, M.-H. Whangbo, andX. G. Gong, Phys. Rev. B , 224429 (2011). H. Xiang, C. Lee, H.-J. Koo, X. Gong, and M.-H.Whangbo, Dalton Trans. , 823 (2013). D. ˇSabani, C. Bacaksiz, and M. V. Miloˇsevi´c, Phys. Rev.B , 014457 (2020). C. Xu, J. Feng, S. Prokhorenko, Y. Nahas, H. Xiang, andL. Bellaiche, Phys. Rev. B , 060404(R) (2020). C. Xu, J. Feng, H. Xiang, and L. Bellaiche, npj Compu-tational Materials , 57 (2018). P. W. Anderson, Phys. Rev. , 2 (1959). S. W. Jang, M. Y. Jeong, H. Yoon, S. Ryee, and M. J.Han, Phys. Rev. Materials , 031001(R) (2019). A. K. Kundu, Y. Liu, C. Petrovic, and T. Valla, ScientificReports , 15602 (2020). A. A. Tsirlin, O. Janson, and H. Rosner, Phys. Rev. B , 144429 (2011). A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, Phys.Rev. B , R5467 (1995). P. Olalde-Velasco, J. Jim´enez-Mier, J. D. Denlinger,Z. Hussain, and W. L. Yang, Phys. Rev. B83