Introduction

The multiferroic materials, exhibiting coupling among different ferroic orders such as ferroelectricity (FE), ferromagnetism (FM), ferroelasticity, and ferrovalleys (FV), etc., have attracted great interest in recent years due to their promising physical properties and potential applications. The spontaneous ferrovalley polarization1 is derived in the 2D FM materials when time-reversal symmetry and space central-inversion symmetry are broken. Two prominent multiferroic materials, i.e. BiFeO3 and TbMnO3 perovskite oxide compounds2,3, have been extensively studied over the past two decades, for their ability to lower energy consumption and simplify device configuration. Recently, the rapid development of two-dimensional (2D) materials, characterized by increased tunability of physical properties compared to bulk systems due to their weak van der Waals (vdW) interlayer interactions4,5,6 has captured attention in both academic and industry communities. Particularly, the investigation of magnetoelectric coupling in 2D materials has been emerging as a promising area in electronic device research7.

However, the practical applications of 2D FE and 2D FM materials have been greatly limited by their scarcity over the past decades. Only a few 2D FE materials, like In2Se38, SnTe9, CuInP2S610, and group-V monolayers (As, Sb, Bi)11,12 have been demonstrated in experiments, leaving 2D multiferroic materials even more sparsely explored. Recently, the concept of sliding-ferroelectrics (SFE) is proposed theoretically, by utilizing the weak interlayer vdW interactions of 2D materials13,14,15. By asymmetrically stacking bilayer 2D materials, a sizable out-of-plane FE polarization can be achieved. This polarization direction can be switched by interlayer sliding between neighboring layers using an external electric field, as confirmed in experimental investigations involving materials such as h-BN16,17, WTe218, β-InSe19, γ-InSe20, 3R-MoS221, among others. This proposal shines a light on the construction of 2D multiferroics with a different coupling mechanism, which have been implemented in a range of materials including 2H-VS222, 1T-FeCl223, VSi2P424, MnBi2Te425, YI226 and MnSe27. The substantial surface-to-volume ratio promotes significant roles for interface engineering and external strain fields in the subtle competition among different ferroic orders in 2D materials. Consequently, external factors like strain28 or electric fields can modulate the coupling among different degrees of freedom—such as electronic charge, spin, and valley—in 2D materials.

In this work, we performed the first-principles calculations on the strain tunable ferroic orders of bilayer R stacking 2D antiferromagnetic VSe2. Our study reveals that a slight out-of-plane pressure induces antiferromagnetism (AFM) to FM interlayer coupling, significantly influencing the coupling between the spin and valley degree of freedom. More interestingly, the splitting of ferrovalley in energy extrema at K+ and K- points of the Brillouin zone shows a linear correlation with interlayer sliding polarization. The strong coupling between the SFE, AFM, and FV allows the mechanical control of ferrovalley as well as the ferromagnetism. This capability holds promise for manipulating exotic physical properties in 2D multiferroics.

Results and discussion

Monolayer properties of VSe2

As shown in Fig. 1, our calculations show that the single layer 2H-VSe2 has a ferromagnetic ground state with a magnetic moment of 1 μB per V atoms. The absence of the imaginary phonon modes in the phonon dispersion indicates its dynamic stability. 2H-VSe2, which belongs to the space group of \({\rm{P}}\bar{6}m2\), intrinsically breaks the spatial inversion (P-) and time (T-) inversion symmetry, resulting in spontaneous ferrovalley polarization1. A valley splitting (ΔE) of 0.11 eV and valley polarization is observed at the K+ and K- points of the Brillouin zone in the calculated band structures (Fig. 1c). By reversing the direction of the magnetic moment, the bandgap of the two energy valleys in the K+ and K- points and valley polarization is reversed accordingly. The magnetocrystalline anisotropy energy (MAE) in the XOY and XOZ planes reveals a preferred magnetization direction in the XOY plane. This coupling between ferrovalley properties and magnetic moment introduces diverse valley physics in R-stacking bilayer antiferromagnetic 2H-VSe2, which lacks both space and time inversion symmetry simultaneously.

Fig. 1: Ferromagnetic-valley coupling in monolayer VSe2.
figure 1

a The ball-and-stick model and the first Brillouin zone of monolayer 2H-VSe2. The dashed lines indicate the unit cell and the shallow yellow dotted surface represents the electric density. Phonon dispersion (b) and spin-polarized band structures (c) of monolayer 2H-VSe2. The solid 3D arrows located across the V atoms denote the direction of magnetic moments. d The magnetocrystalline anisotropy energy (MAE) in XOY and XOZ planes, respectively, where the direction [100] is set as the reference axis.

Full size image

Effect of stackings on ferroelectric properties

As depicted in Fig. 2a–c, there are three different stacking orders for the bilayer VSe2: AA stacking, AB stacking, and BA stacking. AA stacking belongs to \({\rm{P}}\bar{6}m2\) (D3H-1) space group, whose upper layer transitional metal V atoms align to those in the lower layer, resulting in its mirror symmetry Mz. The AB or BA stacking involves a lateral shift of \(\left[\frac{1}{3},\frac{2}{3},0\right]\) or from AA stacking, thus breaking the Mz symmetry. Bilayer AB stacking VSe2 belongs to P3m1 (C3V-1) space group, with an optimized in-plane lattice constant of 3.32 Å and interlayer vdW gap of 3.12 Å. A sizable out-of-plane dipole moment (\({\rm{D}}_{\rm{z}}\)) is formed when the interlayer VSe2 is sliding a short distance from the mirror position between adjacent vdW layers.

Fig. 2: The ferroelectric properties of R-stacking VSe2.
figure 2

The atomic structures (a–c) for AA, AB, and BA stackings of bilayer 2H-VSe2, respectively. The orange and mint green filling colors represent charge accumulation and depletion respectively. Their planar averaged charge differential density (CDD) along z-direction is displayed in (d–f). The red solid arrows represent the directions of electric polarizations in different stackings.

Full size image

Figure 2d–f displays the charge density difference Δρ between the bilayer system and its components for AA, AB, and BA stackings. The plane-average charge density difference in AA stacking reveals a symmetrically distributed charge density, resulting in a zero net interlayer dipole moment. While in AB stacking, a net downward dipole moment (\({{\bf{D}}}_{{\boldsymbol{\downarrow }}}\)) of 1.2 × 10−3 C·m−2 is developed due to the different environment provided by Se atoms around the vdW gap, causing asymmetric distribution of accumulated charge in the out-of-plane direction. A similar analysis is applied to the upward dipole moment (\({{\bf{D}}}_{{\boldsymbol{\uparrow }}}\)) in BA stacking. The dipole moments can be easily switched under external electric fields as the structure of AB stacking can slide to one of its three adjacent BA stackings in space under a small field, giving rise to the vdW sliding ferroelectrics.

Figure 3 shows the stacking energy landscape in the unit cell as a function of the lateral shift. The stacking energy is defined as the energy difference between the shifted configuration in a particular stacking order and the AB stacking order. The lateral shift along the [100] direction corresponds to the sliding of one VSe2 layer along the V-Se bonding direction, while the shift along the \([1\bar{1}0]\) direction corresponds to the nearest V-V direction. We note that AB stacking has the lowest stacking energy and degenerates with BA stacking in energy.

Fig. 3: Effect of interlayer sliding on energy and magnetic ground state.
figure 3

a The energy landscapes of bilayer 2H-VSe2 while moving the upper layer with respect to the fixed lower layer in real space under the paramagnetic state. b The energy difference between interlayer AFM and FM coupling. c, d The energetic profiles along \([1\bar{1}0]\) and [100] directions (denoted as orange and blue dashed lines in (a, b)) for (a, b), respectively.

Full size image

Employing the nudged-elastic band (NEB) method, we explored the transition path along the \([1\bar{1}0]\) direction. Along this path, an energetically lower transition state than AA (calculated as 134 meV), labeled as ST, is observed. ST corresponds to a fractional lateral shift of \(\left[\frac{1}{2},\frac{1}{2},0\right]\) with respect to the AA stacking. As shown in Fig. 3c, to switch the dipole moment from the \({{\bf{D}}}_{{\boldsymbol{\downarrow }}}\) to \({{\bf{D}}}_{{\boldsymbol{\uparrow }}}\) by lateral shifting, it needs to overcome a small energy barrier of 9.9 meV per unit cell across the ST point. This barrier is lower compared to the 2D ferroelectric In2Se3 (about 66 meV per unit cell)29 and traditional 3D perovskite PbTiO3 (about 240 meV per unit cell)30, but comparable to h-BN (about 4 meV per unit cell)13. This suggests that the polarization in AB stacking bilayer VSe2 is stable in the ground magnetic state, and the small sliding barrier can be easily overcome by an external electric field, as demonstrated in bilayer h-BN SFE by Wu13.

We also calculated the interlayer exchange energy as a function of the lateral shift in Fig. 3b. The interlayer exchange energy is defined as the energy difference between the FM states and interlayer AFM states in bilayer VSe2. Interestingly, the bilayer VSe2 prefers the AFM ground state, showing a comparatively lower interlayer exchange energy of −0.82 meV in the AB (BA) stacking. The interlayer exchange energy is much larger in the AA and ST, measuring −8.35 meV and −1.36 meV, respectively. The ultralow exchange energy in the AB stacking also suggests that it could be easily tuned by external field.

Effect of stackings on ferrovalley properties

Combining sliding ferroeletric and interlayer exchange coupling, it could easily create four multiferroic states in bilayer VSe2, i.e. \({{\bf{D}}}_{{\boldsymbol{\uparrow }}}{{\bf{M}}}_{{\boldsymbol{\uparrow }}{\boldsymbol{\downarrow }}}\), \({{\bf{D}}}_{{\boldsymbol{\downarrow }}}{{\bf{M}}}_{{\boldsymbol{\uparrow }}{\boldsymbol{\downarrow }}}\), \({{\bf{D}}}_{{\boldsymbol{\uparrow }}}{{\bf{M}}}_{{\boldsymbol{\downarrow }}{\boldsymbol{\uparrow }}}\) and \({{\bf{D}}}_{{\boldsymbol{\downarrow }}}{{\bf{M}}}_{{\boldsymbol{\downarrow }}{\boldsymbol{\uparrow }}}\), where \({{\bf{M}}}_{{\boldsymbol{\downarrow }}{\boldsymbol{\uparrow }}}\) denotes the magnetic moment of the top and bottom VSe2 layer pointing downward and upward. The band structures of \({{\bf{D}}}_{{\boldsymbol{\uparrow }}}{{\bf{M}}}_{{\boldsymbol{\uparrow }}{\boldsymbol{\downarrow }}}\) and \({{\bf{D}}}_{{\boldsymbol{\downarrow }}}{{\bf{M}}}_{{\boldsymbol{\downarrow }}{\boldsymbol{\uparrow }}}\) appear indistinguishable as shown in Fig. 4a and e, so is \({{\bf{D}}}_{{\boldsymbol{\downarrow }}}{{\bf{M}}}_{{\boldsymbol{\uparrow }}{\boldsymbol{\downarrow }}}\) and \({{\bf{D}}}_{{\boldsymbol{\uparrow }}}{{\bf{M}}}_{{\boldsymbol{\downarrow }}{\boldsymbol{\uparrow }}}\). For a particular state featuring both SFE and AFM order, the energy valleys at the adjacent hexagonal Brillouin zone K+ and K- no longer remain degenerate due to the breaking of time reversal symmetry and inversion symmetry. The calculated Berry curvatures shown in Fig. 4c and f have opposite signs in value for the K+ and K- valleys, suggesting that there are different pseudo-magnetic fields in these valleys. Thus the electrons in adjacent K valley will have an opposite magnetic force and cause the different spin current flow in the top and bottom surface due to the anomalous valley Hall effect. By comparing the electronic structures of different SFE and AFM states in Fig. 4, we found that the band structure is related to Mirror symmetry Mz breaking during sliding, indicating an interconnected relationship between valley polarization and ferroelectric polarization.

Fig. 4: Valley polarization properties under different multiferroic couplings.
figure 4

The spin-projected band structures (left panel) and its schematic graph (right panel) of VSe2 for AB stacking under a head-to-head (a)/tail-to-tail (b) interlayer magnetic moment. Panels d, e is the same conditions for BA stacking. The corresponding calculated Berry curvature for (c) AB stacking and (f) BA stacking in head-to-head AFM configuration. The short dashed and solid lines in the schematic graph represent the contributions from the bottom and upper layer, respectively. The spin-up and spin-down band structures are indicated by red and blue color.

Full size image

There is a strong coupling between spin and valley degree of freedom in the bilayer 2H-VSe2, as evident in Supplementary Fig. 1, where the V atoms in different layers contribute different energy extrema in the valence band maximum (VBM) and conduction band minimum (CBM). For instance, considering the \({{\bf{D}}}_{{\boldsymbol{\downarrow }}}{{\bf{M}}}_{{\boldsymbol{\downarrow }}{\boldsymbol{\uparrow }}}\) configuration in Fig. 4a, the upper VSe2 layer contributes the blue spin-down projected band in the CBM while the bottom layer provides the red spin-up projected VBM for AB stacking (also shown in layer resolved band structures in Supplementary Fig. 1a), confirming the layer-dependent spin polarization24,31,32. While in the \({\rm{K}}_{\pm }\) valleys of BA stacking with the \({{\bf{D}}}_{{\boldsymbol{\uparrow }}}{{\bf{M}}}_{{\boldsymbol{\downarrow }}{\boldsymbol{\uparrow }}}\) configuration (Fig. 4c and Supplementary Fig. 1b), the spin-down VBM and spin-up CBM characteristic bands are contributed by upper and bottom layers, respectively. Conversely, in AA stacking, as plotted in Fig. 5c, \({\rm{K}}_{\pm }\) valley remain degenerate, exhibiting zero energy difference in the VBM due to the protection of the joint symmetry \(\hat{{\rm{O}}}\equiv \hat{{{\rm{M}}}_{{\rm{z}}}}\hat{{\rm{T}}}\) in space. Furthermore, when the bilayer VSe2 is in the FM state, the spin-polarized bands do not have layer dependence in the FM state, and the VBM and CBM are mainly contributed from the red spin-up and blue spin-down bands, as shown in Supplementary Fig. 1d–f.

Fig. 5: Multiferroic properties under different operations.
figure 5

a Electric polarization D variation as a function of interlayer displacement. Red solid line circles the special points mentioned in our context. b Valley split value (UE) of bilayer 2H-VSe2 plotted against ferroelectric polarization D and its linear fitted line. The fitting slope is 73.13 eV·m2·C−1. c The transition barrier for AB-BA switching along the direction \([1\bar{1}0]\) for VSe2. d Energy difference ΔE2 = EAFM-EFM between interlayer AFM and FM coupling with respect to the c-axis strain (defined as \(\frac{c}{{c}_{0}-1}\)).

Full size image

How do the stacking orders and magnetic configurations affect the valley degree of freedom? To answer this question, we enlarged the spin-projected electronic band structures of the \({{\bf{D}}}_{{\boldsymbol{\uparrow }}}{{\bf{M}}}_{{\boldsymbol{\downarrow }}{\boldsymbol{\uparrow }}}\) \({\rm{and}}\) \({{\bf{D}}}_{{\boldsymbol{\downarrow }}}{{\bf{M}}}_{{\boldsymbol{\downarrow }}{\boldsymbol{\uparrow }}}\) in Supplementary Fig. 2a and b, where the \({{\bf{M}}}_{{\boldsymbol{\downarrow }}{\boldsymbol{\uparrow }}}\) AFM bilayer VSe2 has a “head to head” magnetic moment in each layer. The \({{\bf{D}}}_{{\boldsymbol{\uparrow }}}{{\bf{M}}}_{{\boldsymbol{\downarrow }}{\boldsymbol{\uparrow }}}\) bilayer VSe2 has a spin splitting of 81.7 meV and 49.9 meV in the K+ and K- valleys for the lowest two CBM, the corresponding values are 6.6 meV and 164.4 meV for the highest two VBM (also shown in Supplementary Fig. 2a, b). While in the \({{\bf{D}}}_{{\boldsymbol{\downarrow }}}{{\bf{M}}}_{{\boldsymbol{\downarrow }}{\boldsymbol{\uparrow }}}\) configuration, the spin splittings have similar values, but in opposite \({{\rm{K}}}_{\pm }\) valleys. However, as shown in Supplementary Fig. 2c, there is no such a spontaneous ferrovalley in the mirror symmetry AA stacking, although it still has the spin dependent \({{\rm{K}}}_{\pm }\) valley degree of freedom. A calculated energy difference (UE) of 79.1 meV for the VBM between the K+ and the K- valley is found in the AB stacking. This energy difference between the VBM positions in K+ and K- valley is comparable with other monolayer magnetic materials, e.g. 105 meV in H-FeCl233, 95 meV in ScI234, and 156 meV in NbS235, but is considerably large compared with other bilayer systems such as bilayer T-FeCl2 (4.7 meV)23 and bilayer YI2 (24 meV)26. The direction of the dipole polarization can influence the energy extrema of VBM at the K valleys, we depicted the corresponding valence band alignment schematically in Fig. 4 for the AB and BA stacking orders. Therefore, the layer-dependent spin-polarized band structure can be readily modulated by sliding ferroelectric orders. Since the valley degree of freedom can be detected by the circular polarized light, the spin polarization, polarized photon, and ferroelectric polarizations are related to each other in the multiferroic bilayer VSe2. Notably when the AB stacking bilayer VSe2 is in the FM states (Supplementary Fig. 2d), the ferrovalley remains with valley splitting of 9.4 meV between the K+ and K- points. However, different from the monolayer FM VSe2, the VBM is located at Γ point with a much smaller bandgap of 0.034 eV (Γ to M) in bilayer FM VSe2.

k·p model analysis

To quantitatively analyze the influence of sliding polarizations on the ferrovalley and spin polarization, we employ an effective k·p model to discuss the interactions of different degrees of freedom in the model by partitioning their effect on the spin-polarized electronic band structures. The effective Hamiltonian of bilayer VSe2 can be expressed as

$${{\rm{H}}}_{{\rm{k}}}={{\rm{I}}}_{2}\otimes \left[\begin{array}{cc}{{\rm{H}}}_{{\rm{k}}}^{{\rm{u}}} & {{\rm{H}}}_{\perp }\\ {{\rm{H}}}_{\perp } & {{\rm{H}}}_{{\rm{k}}}^{{\rm{b}}}\end{array}\right]$$
(1)

Here, I2 is a 2 × 2 identity matrix. \({{\rm{H}}}_{\perp }\) represents the interlayer hopping term, which is ignored in our analysis due to the weak vdW interaction. \({{\rm{H}}}_{{\rm{k}}}^{{\rm{u}}({\rm{b}})}\) represents the Hamiltonian of the upper (bottom) layer. This model comprises the nearest hopping term, spin-orbit coupling (SOC), exchange interaction, and electric polarization arising from interlayer sliding:

$${{\rm{H}}}_{{\rm{k}}}^{{\rm{u}}}=\left[\begin{array}{cc}\frac{\Delta }{2}+\varepsilon +{\rm{\tau }}{\rm{s}}{\lambda }_{{\rm{c}}}-{\rm{s}}{m}_{{\rm{c}}}+\frac{{\rm{U}}_{{\rm{E}}}}{2} & {{\rm{t}}}_{12}\left({\rm{\tau }}{q}_{x}-{\rm{i}}{q}_{y}\right)\\ {{\rm{t}}}_{12}\left({\rm{\tau }}{q}_{x}+{\rm{i}}{q}_{y}\right) & -\frac{\Delta }{2}+{\rm{\varepsilon }}+{\rm{\tau }}{\rm{s}}{\lambda }_{{\rm{v}}}-{\rm{s}}{m}_{{\rm{v}}}+\frac{\rm{U}_{{\rm{E}}}}{2}\end{array}\right]$$
(2)
$${{\rm{H}}}_{{\rm{k}}}^{{\rm{b}}}=\left[\begin{array}{cc}\frac{\Delta }{2}+\varepsilon +{\rm{\tau }}{\rm{s}}{\lambda }_{{\rm{c}}}+{\rm{s}}{m}_{{\rm{c}}}-\frac{{\rm{U}}_{{\rm{E}}}}{2} & {{\rm{t}}}_{12}\left({\rm{\tau }}{q}_{x}-{\rm{i}}{q}_{y}\right)\\ {{\rm{t}}}_{12}\left({\rm{\tau }}{q}_{x}+{\rm{i}}{q}_{y}\right) & -\frac{\Delta }{2}+\varepsilon +{\rm{\tau }}{\rm{s}}{\lambda }_{{\rm{v}}}+{\rm{s}}{m}_{{\rm{v}}}-\frac{{\rm{U}}_{{\rm{E}}}}{2}\end{array}\right]$$
(3)

Where ∆ represents the band gap at K+ (K-) valley, ε is the energy correction related to the fermi level, \({\rm{\tau }}=\pm 1\) is the valley index, t is an effective nearest-neighbor hopping integral, and q = k- K denotes the relative momentum vector with respect to the K+ (K-) valley. \({\rm{s}}=\pm 1\) is the spin index where its sign represents spin-up and spin-down respectively. \({\lambda }_{{\rm{c}}({\rm{v}})}={\rm{E}}_{{\rm{c}}({\rm{v}})\uparrow }-{\rm{E}}_{{\rm{c}}({\rm{v}})\downarrow }\) is the spin splitting at CBM (VBM) in the single-layer TMDs derived from SOC effect. \({m}_{{\rm{c}}({\rm{v}})}={\rm{E}}_{{\rm{c}}({\rm{v}})\downarrow }-{\rm{E}}_{{\rm{c}}({\rm{v}})\uparrow }\) is the spin splitting at CBM (VBM) in the single layer TMDs due to exchange interaction, and \({\rm{U}}_{{\rm{E}}}\) is the induced dipole energy of each layer derived from the interlayer charge transfer. For the \({\bf {D}}_{\uparrow }{\bf{M}}_{\downarrow \uparrow}\) configuration, the corresponding eigenvalues are obtained by performing diagonalization operations on the Hamiltonian:

$$\begin{array}{c}{\rm{E}}\left({\psi }_{{\rm{u}}}^{\tau }\right)=\frac{1}{2}\left[2\varepsilon -{\rm{\tau }}\left({\lambda }_{{\rm{c}}}+{\lambda }_{{\rm{v}}}\right)+\left({m}_{{\rm{c}}}+{m}_{{\rm{v}}}\right)\pm \sqrt{{\left({{\Delta}^{\prime} }_{{\rm{u}}}\right)}^{2}+4{{\rm{t}}}_{12}^{2}\left({{\rm{\tau }}}^{2}{q}_{x}^{2}+{q}_{y}^{2}\right)}\right]+\frac{{\rm{U}}_{{\rm{E}}}}{2}\\ {\rm{E}}\left({\psi }_{{\rm{b}}}^{\tau }\right)=\frac{1}{2}\left[2\varepsilon +{\rm{\tau }}\left({\lambda }_{{\rm{c}}}+{\lambda }_{{\rm{v}}}\right)+\left({m}_{{\rm{c}}}+{m}_{{\rm{v}}}\right)\pm \sqrt{{\left({{\Delta}^{\prime} }_{{\rm{b}}}\right)}^{2}+4{{\rm{t}}}_{12}^{2}\left({{\rm{\tau }}}^{2}{q}_{x}^{2}+{q}_{y}^{2}\right)}\right]-\frac{{\rm{U}}_{{\rm{E}}}}{2}\end{array}$$
(4)

Here, the gap parameter \({{\Delta}^{\prime} }_{{\rm{u}}\left({\rm{b}}\right)}=\Delta -\left(+\right){\rm{\tau }}\left({\lambda }_{{\rm{c}}}-{\lambda }_{{\rm{v}}}\right)+\left({m}_{{\rm{c}}}-{m}_{{\rm{v}}}\right)\) reflects the band gap of each sub-band contributed by different layers. We noted that the interlayer sliding can easily modulate the relative energy difference between the two sub-bands. The difference between VBM of the two sub-bands is characterized by \({\rm{U}}_{{\rm{E}}}\), expressed as \({\rm{E}}\left({\psi }_{{\rm{u}}}^{1}\right)-{\rm{E}}\left({\psi }_{{\rm{b}}}^{-1}\right)={\rm{U}}_{{\rm{E}}}\). This indicates that the relative altitude of VBM for spin-up and spin-down can be singly tuned by the out of plane (OOP) polarization associated with the interfacial sliding. Since the CBM of K+ and K- valleys have similar values, so the \({\rm{U}}_{{\rm{E}}}\) quantity could be considered as a rough estimation of bandgap difference between K+ and K- valleys. It’s noteworthy that the OOP polarization does not affect the shapes of the band structures. Referring to the band alignment in Fig. 4 (also Supplementary Fig. 2), we can define the direction of valley polarization based on the sign of \({\rm{U}}_{{\rm{E}}}\). The switching from positive to negative valley polarization can be realized by interlayer sliding along the direction \(\left[\frac{2}{3},-\frac{2}{3},0\right]\) or \(\left[-\frac{2}{3},\frac{2}{3},0\right]\). By comparing the band structures in \({{\bf{D}}}_{{\boldsymbol{\uparrow }}}{{\bf{M}}}_{{\boldsymbol{\downarrow }}{\boldsymbol{\uparrow }}}\), \({{\bf{D}}}_{{\boldsymbol{\downarrow }}}{{\bf{M}}}_{{\boldsymbol{\downarrow }}{\boldsymbol{\uparrow }}}\), \({{\bf{D}}}_{{\boldsymbol{\uparrow }}}{{\bf{M}}}_{{\boldsymbol{\uparrow }}{\boldsymbol{\downarrow }}}\) and \({{\bf{D}}}_{{\boldsymbol{\downarrow }}}{{\bf{M}}}_{{\boldsymbol{\uparrow }}{\boldsymbol{\downarrow }}}\) states in Fig. 4, it becomes evident that external fields such as shear sliding, electric field, or magnetic field can be exerted to realize a “four-state switching” between different arrangements of valley polarization and SFE polarization, rendering a ferroelectricity-valley coupling.

Different means of modulations

The interlayer sliding plays a dual role in directing the electric and valley polarizations simultaneously. Our k·p model reveals a link between the difference in valley VBM at K+ and K- and the interlayer electric polarization. In Fig. 5a, the induced electric dipole polarization is displayed as a function of interlayer displacement during the sliding ferroelectric transition. The electric polarization shows a centrosymmetric tendency with a zero polarization point located at the saddle point ST, whose transition barrier and atomic structure are shown in Fig. 5c. Interestingly, we also calculated the energy difference of the VBM at K+ and K- valleys as a function of the sliding displacement, which shows a similar tendency to the electric polarization. Figure 5c plots the scatter diagram of the energy difference value between K+ and K- valleys UE with respect to the electric polarization. The fitted function shows that the UE is linearly related to the electric polarization with a slope factor of 73.13 eV·m2·C−1. This indicates a strong coupling between interlayer electric polarizations and valley polarizations in the SFE bilayer AFM 2H-VSe2.Furthermore, as depicted in Fig. 3d, there is a giant magneto-electric coupling in the AB stacking bilayer 2H-VSe2 due to the small interlayer exchange energy. The magnetism can be easily controlled by external fields such as strain or magnetic field. The plot in Fig. 5d demonstrates the interlayer strain-dependent magnetic orders for 2H-VSe2, where the interlayer AFM ground state switches to the FM ground states under a 5.07% compressive strain. Based on the calculated elastic constants of C33 in Table 1, a small pressure of 3.2 Gpa is sufficient to induce the AFM to FM transition in bilayer 2H-VSe2, leading to completely different valley responses shown in Supplementary Fig. 2. It is also noted that there is a semiconductor to metal transition when the compressive strain is gradually increased in the bilayer 2H-VSe2, as shown in Supplementary Fig. 5.

Table 1 Elastic constants
Full size table

Other candidate materials and outlooks

Our approach can extend to other 2D VX2 (X = S, Te) magnetic bilayer materials. By replacing the Se element with Te and S element, we show that VTe2 and VS2 also exist similar linear ferroelectricity-valley coupling as observed in VSe2, suggesting the potential for simultaneous SFE and valley polarization switching. As shown in Fig. 6a and c, the switching barriers during the same sliding path for VS2 and VTe2 are 4.88 and 5.97 meV, respectively. Figure 6b and c exhibit the scatter diagrams for the valence band energy difference (UE) between K+ and K- valley as a function of the electric polarizations D. The ratio α of UE to D shows similar values across VSe2, VS2, and VTe2. Their optimized structure parameters, ferroelectric polarization, valley splitting and UE are listed in Table 2. There exist slight differences in polarization values due to varying charge redistributions among systems. Besides, the electronic band structures under different interlayer AFM configurations ensure the “four state switching”22,26, as shown in the Supplementary Fig. 3 and 4. The strong linear correlation between UE and electric polarization also provides an alternative method to measure tiny SFE polarization by chiral circularly polarized light. Similar to VSe2, the bilayer VTe2 also has an AFM ground state, and a compressive strain of 1.8% can cause the transition from AFM to FM, as shown in Fig. 6f. In contrast, the bilayer VS2 has the FM ground state (Fig. 6c). Based on the calculated elastic constants for VS2 and VTe2 in Table 1, predictions suggest the critical strain for the AFM to FM transition in bilayer VS2 and VTe2 to be −2.02 GPa (where a negative sign implies tensile strain) and 1.20 GPa, respectively. Based on the above strain-induced AFM to FM phase transition, a multiple external field tunable spin valve device is proposed in Fig. 7. The bilayer AFM 2H-VSe2 is placed inside a transparent diamond anvil, which could allow the circularly polarized light to penetrate through the anvil. There are both spin-up and spin-down electrons excited under different polarized light when there is no strain, and the flow of spin electrons results in the different spin currents on the left and right edges due the anomalous valley Hall effect. Once the AFM to FM phase transition is induced by an external strain, only one kind of spin current is generated by the polarized light as both K+ and K- valleys generate the same spin electrons. So it could be used to generate different spin currents by the strain engineering. Additionally, applying gate voltage can also induce valley polarization switching, as shown in Fig. 4a and d, which can also be used to switch the flow direction of spin currents.

Fig. 6: Multiferroic properties in VS2 and VTe2.
figure 6

The energy barrier of electric polarization switching for VS2 (a) and VTe2 (d) in the direction \([1\bar{1}0]\). b, e Corresponding scatter diagrams and its linear fitting for UE-D relationship. Energy difference ΔE2 = EAFM-EFM between interlayer AFM and FM coupling with respect to the c-axis strain for VS2 (c) and VTe2 (f).

Full size image
Table 2 Parameters for multiferroic coupling
Full size table
Fig. 7: Strain-tunable spin valve device.
figure 7

The Spin valve devices are designed by manipulating strain and gate voltage to control valley polarization and the flow of spin current.

Full size image

In conclusion, based on first-principles density functional theory calculation and the effective k·p model, we unveiled the intrinsic coupling between the ferrovalley and the sliding ferroelectric polarizations in bilayer AFM 2H-VSe2. The bilayer AFM VSe2 exhibits a layer-dependent spin-polarized band structure, resulting in an energy difference between the VBM in K+ and K- valleys. The breaking of mirror symmetry plays a critical role in the valley-electric polarization coupling, where dipole polarization and valley polarization arise spontaneously from interlayer sliding. Our DFT calculations highlight a strong linear relationship between the energy difference of VBM among different K valleys and electric polarization, offering a means to measure electric polarization. Similar phenomena were observed in VX2 (X = S and Te) compounds. The robust coupling between strain, ferromagnetic orders, electric, and valley polarization enables the manipulation of functional properties in bilayer VX2 systems through strain engineering, promising the design of multifunctional devices.

Methods

Electronic structure calculations

The first-principles density functional theory (DFT) calculations were carried out in the projector augmented wave method36 via the Vienna ab initio simulation package (VASP)37,38. The exchange-correlation function was described by the generalized gradient approximation (GGA) within the Perdew-Burke-Ernzerhof (PBE) formalism39. Structures were considered as fully relaxed when the maximum force between atoms is less than 0.01 eV Å−1 and an energy cutoff of 500 eV was used for the plane wave expansion. The first Brillouin zone was sampled by a Γ-centered Monkhorst–Pack k-point mesh of 15 × 15 × 1. The strong correlation effect of the d electron of V atoms was considered via the on-site Hubbard correction PBE + U method40 with the effective U = 1.0 eV41 in the structural optimization and electronic self-consistent calculation. The vacuum space of 20 Å was introduced to avoid interactions between the adjacent slabs. DFT-D3 method was used to consider the van der Waals correction in all of our bilayer VX2. The phonon dispersion over a 4 × 4 × 1 supercells were calculated in VASP through PHONOPY package42,43 using harmonic approximation.