Introduction

Multiferroic materials are a special class of multifunctional materials that possess broken space inversion symmetry, switchable polarization under external fields and magnetism. They have been reported to display varieties of intriguing physical phenomena including, but not limited to, large piezoelectric coefficients, dielectric tunability, complex structural transitions, nontrivial dipole and spin textures and magnetoelectric coupling effects, which lead to broad applications such as actuators, sensors, field-effect transistors, tunneling junctions and emergent electronics1,2,3,4,5,6. For applications in miniaturized memory or logic devices, it is highly desirable to have a multiferroic material as thin as a few unit-cell thickness with sharp and clean surface, robust out-of-plane (OP) electric polarization, magnetic orderings and strong magnetoelectric coupling sustainable above room temperature (RT). In this regard, multiferroicity in 2D van der Waals (vdW) materials consisting of a single or few layers offers promising alternatives. 2D materials inspired by single-atomic-layer graphene and BN has flourished and drawn great attention in the past decade7. Unfortunately, only a few 2D materials have been experimentally verified to be intrinsic ferroelectrics, and 2D intrinsic ferroelectric materials with OP polarization are very scare so far. For instance, in-plane (IP) ferroelectricity is observed in monolayer SnTe8,9,10, while OP polarization has been reported in few-layer CuInP2S6 nanoflakes11, In2Se3 with thickness as thin as three layers12,13, MoTe2 down to the single atomic layer14,15 and Bi monolayer16. There are also some other computational predictions on 2D intrinsic ferroelectrics which still await experimental confirmation17,18,19. On the other hand, 2D intrinsic multiferroic materials with strong magnetoelectric effects are even more scare. CuCrP2S6, Ti3C2Tx MXene and NiI2 have only been recently confirmed to show both ferroelectricity and magnetism20,21,22. In addition, some approaches have been proposed to induce multiferroicity in 2D materials23,24,25,26,27,28,29,30,31. Among these, polar stacking of 2D vdW materials is a straightforward and efficient way32,33,34. However, stacking-based 2D multiferroics typically display a notably small polarization and minimal intrinsic coupling between ferroelectricity and magnetism due to their inherent separation in origins32,33,34,35,36,37,38,39. Then the question remains: is it possible to achieve a large electric polarization (compared to known stacking 2D ferroelectrics) and strong magnetoelectric coupling simultaneously in a single-phase vdW material under polar stacking?

In this work, using first-principles calculations, we show that bilayer Hf2S presents a much larger OP electric polarization by interlayer stacking despite of its metallicity, as compared to most stacking 2D ferroelectrics. This large electric dipole is revealed to stem from the delocalized interlayer electrons. In addition, there is an unexpected emergence of net magnetization which can be ascribed to Hf atoms at the surface layers, given the fact that the bulk phase is nonmagnetic while the monolayer phase is antiferromagnetic. More importantly, a strong magnetoelectric coupling is identified that ensures the reversal of magnetization with the switching of electric dipoles. The coexistence of polar electric dipole, magnetization and metallicity in one 2D material makes it a rare multiferroic metal. Our work presents an alternative way of searching for multiferroics with OP polarization and strong magnetoelectric couplings in ultrathin materials.

Results and discussion

Structural and ferroelectric properties of Hf2S bilayer

Hf2S is a known example of layered electride materials which are characterized by the presence of electrons trapped within the vdW gap region between layers (also referred as interlayer anionic electrons)40. Bulk Hf2S can be viewed as consisting of two layers of which one is rotated by 180° with respect to the other one. The upper panel of Fig. 1a gives the bulk structure (with space group \(P{6}_{3}/{mmc}\)). For each monolayer, each S atom is sandwiched by the top and bottom Hf layers and bonds to six Hf atoms. Neither the bulk nor the bilayer with the same stacking order as in the bulk give a net electric dipole, due to the presence of spatial inversion symmetry. However, with the polar stacking shown in Fig.1b, bilayer Hf2S permits a net electric dipole. In this polar stacking, the nominal anion S atoms from the lower layer are situated right below the cations Hf from the upper layer, while the cations Hf from the lower layer are located right below the hollow centers from the upper layer (see the top view in Fig.1b). Such polar stacking gives a large electric dipole (0.014 \({\rm{e}}{\text{\AA }}\)/u.c.) pointing downwards, which is much larger than that (0.003 \({\rm{e}}{\text{\AA }}\)/u.c.) of the polar stacking bilayer BN37. Phonon dispersion in Supplementary Fig. 1 indicates its dynamical stability.

Fig. 1: Cystal structures and ferroelectricity of bulk and the polar stacking bilayer of Hf2S.
figure 1

a The crystal structure for bulk Hf2S: side view (left) and top view (right). b The polar stacking of bilayer: side view (left) and top view (right). c Isosurface for differential charge density between the polar stacking bilayer Hf2S and two single monolayers with the same atomic positions. The value is set to 20% of the maximum. The markers X2, X1d and X1u points present the interlayer anionic electron charge centers in bilayer Hf2S. d The homogeneous switching barrier (left axis) and the corresponding polarization (right axis) from NEB calculation.

Full size image

To understand how the electric dipole forms, the charge redistribution is shown in Fig. 1c under such polar stacking, by plotting the differential charge density between the bilayer and two monolayers (see caption for details). Charge accumulation (shown in yellow color) occurs at X2 along the Hf-S vertical segment (denoted by dashed line) and around the interlayer anionic center X1u, while charge depletion occurs around the interlayer anionic center X1d. Notably, the charges transferred from around X1d to X2 lead to the generation of a downward electric dipole and the disappearance of the interlayer anionic electron at X1d (see Supplementary Fig. 2 for comparison). In the case of the polar stacking of other 2D bilayer materials reported in the literature such as BN, the electric dipole comes from the distorted electron cloud along the cation-anion vertical segment (the one denoted by the dashed line)37. Here, the large electric dipole in bilayer Hf2S additionally benefits from the interlayer anionic electron which is delocalized and able to relocate during interlayer stacking shifting (see Supplementary Fig. 2), which may not be the case for other layered electrides with strongly localized interlayer anionic electron (e.g., see Gd2C in Supplementary Fig. 5). To demonstrate that it can be ferroelectric, that is the electric dipole can reverse under an external electric field, an upper bound for the switching barrier is then estimated by means of the Nudged Elastic band (NEB) method41, which is 125 meV/f.u., as shown in Fig. 1d. This implies that the electric dipole may be reversible through interlayer sliding from the perspective of the small switching barrier (see also Supplementary Table 3).

Interestingly, Fig. 2a, b provide the atom-resolved band structure and density of states for the polar stacking Hf2S bilayer, which clearly shows that the system remains metallic. Contrary to the common belief that ferroelectricity and metallicity are mutually exclusive, we show here that bilayer Hf2S is a potential “ferroelectric” metal with coexisting electric dipoles and conducting free carriers. By further inspecting the band structure and density of states, one can see that there are more electronic states contributed from the top and down surface regions than the regions near the vdW gap at the Fermi level, which means that a large portion of free carriers is located at surfaces. While from Fig. 1c, it has been already demonstrated that the ferroelectric dipoles mostly come from the electron redistribution near the vdW gap. The space separation of free carriers and ferroelectric dipoles may lead to weak screening, which sustains their coexistence42.

Fig. 2: Electronic properties of polar stacking bilayer Hf2S.
figure 2

a The projected spin-polarized electronic band structure. b The projected density of states for Hf1 (dotted lines) and Hf2 (solid lines). c Schematic of Hf1 5d orbitals’ splitting due to crystal field, band hybridization and exchange interaction. d The projected density of states of 5d orbitals from Hf1 (black dotted line) and Hf2 (red solid line) in nonmagnetic case.

Full size image

Electronic and magnetic properties of Hf2S bilayer

Now, let us take a closer look at the electronic properties of the polar stacking bilayer Hf2S. As compared to the nonmagnetic and antiferromagnetic nature of the bulk and nonpolar stacking bilayer Hf2S (see Supplementary Fig. 3), both the band structure and density of states indicate that the polar bilayer Hf2S is a system with a net magnetization instead. As similar to the monolayer case, the local magnetic moment of Hf is thought to be of Stoner type and to result from the large 5d orbital-derived states of the surface Hf1,2 atoms (denoted by Hf1 on the bottom surface and Hf2 on the top surface) around the Fermi level (see Supplementary Fig. 4), which is in-line with previous reports43. Due to the fact that Hf 5d orbital is quite delocalized (see Fig. 2a,b), the on-site Coulomb interactions between them and their exchange splittings are likely to be not large. As a result, no band gap opens at the Fermi level. The calculated magnetic moments of the surface Hf atom are around 0.43 μB and 0.08 μB for Hf1 and Hf2, respectively, with their directions remaining antiparallel. The net magnetism in this polar bilayer Hf2S is found to be around 0.35 μB/f.u., which mainly comes from the inequivalence of Hf1 and Hf2 atoms at the two opposite surfaces. Figure 3a shows the isosurface of the magnetization density for the polar stacking bilayer. The magnetization density from up and down spin channels are noticeably inequivalent and located around Hf atoms at the two opposite surfaces, as aforementioned. To better understand the difference in magnitude between the magnetic moments of Hf on the two opposite surfaces, a detailed analysis on the decomposition of the electronic structures is warranted. The atom-resolved band structures and density of states indicate that the Hf 5d orbital-derived states are split into three subsets a1g (dz2), e1g (dx2-y2, dxy) and e2g (dxz, dyz) with the trigonal crystal field, as shown in Fig. 2c. Note that the Hf 5d orbitals are not localized but rather quite dispersive. For instance, the Hf 5dxy, 5dxz-derived states extend almost all over the energy window (over 8 eV in Fig. 2a,b). Therefore, the three subsets will further hybridize and split, as shown in the schematic diagram in Fig. 2c. The occupations among these states are thus different for Hf1 and Hf2 atoms. Due to the electron redistribution occurring mainly in the vdW gap in the polar stacking bilayer Hf2S, the generation of ferroelectric dipoles then leads to a potential difference between Hf1 and Hf2 atoms. For the dipole moment pointing downward, the energy levels derived from the orbitals (particularly d orbitals) of Hf2 atoms will shift higher. Figure 2d gives the projected density of states from Hf1 and Hf2 atoms on the surfaces for the nonmagnetic case, which indicates that there is a decrease in the density of states of Hf2 at the Fermi level, as compared to that of Hf1 (see also Fig. S4b). As mentioned above, the magnetism of Hf atom is deemed to be of Stoner type. Such decrease of density of states suggests a weakening of the Stoner magnetism, which leads to a further decrease of the exchange splitting and occupation redistribution among Hf2 d-derived orbitals. Figure 2b clearly shows that there is more (less) occupation of e1g (e2g) subset in the spin up (down) channel for Hf2, which ultimately leads to a reduction of the spin magnetic moment of Hf2 and a sizeable net spin magnetic moment that makes the system ferrimagnetic (see Fig. 3a). Such delicate balance between band filling and the Stoner magnetism instability under the influence of the inherent electric field from the polar stacking is proposed to underlie the strong magnetoelectric effect discussed later.

Fig. 3: Magnetic properties of Hf2S.
figure 3

a Isosurface of the magnetic charge density. The value is set to 20% of the maximum. b Simulated magnetization versus temperature.

Full size image

To further examine the stability of the ferrimagnetism present in the system, conventional energy mapping method is adopted to extract the exchange parameters between neighboring magnetic Hf atoms (see Supplementary Fig. 6 for more details). Considering that the magnetic moments on the Hf3,4 (near the vdW gap) and Hf2 ions are comparably small and negligible, magnetic interactions for as far as third nearest neighbors on the bottom surface are tested to be sufficient. With the extracted exchange parameters, a Curie temperature of around 248 K is estimated by the Monte Carlo simulations (Fig. 3b and see Supplementary Fig. 6 for details).

Magnetoelectric coupling and its mechanism

As shown above, the polar stacking which leads to ‘ferroelectricity” also generates net magnetization in bilayer Hf2S, which suggests a possible correlation between ferroelectricity and ferrimagnetism in this material, and thus makes the polar stacking Hf2S bilayer a potential multiferroic-type material with magnetoelectric coupling. To demonstrate it, we first turn to the evolution of magnetism during the switching of ferroelectric polarization. As it has been shown in Fig. 1d, the switching through interlayer sliding is possible in the studied material. Figure 4a shows a schematic of the proposed switching process of ferroelectric dipoles and the accompanied evolution of magnetism. At stage I, the system is in one of its degenerate ground states, i.e., with electric polarization pointing downward and magnetization pointing to the right. At stage II, with the reduction of the electric polarization by sliding the top layer to the left, the magnetic moment of Hf2 increases, while that of Hf1 remains almost the same. Consequently, the net magnetic moment decreases. At stage III, with the increase of the magnetic moment of Hf2 and Hf4, the overall net magnetization further decreases and then changes its direction and points to left. At stage IV, the ferroelectric polarization disappears, and the magnetic moment of Hf2 increases to be the same as that of Hf1 in magnitude but antiparallel in direction, and therefore no net magnetization exists.

Fig. 4: Magnetoelectric coupling.
figure 4

a Snapshots of electric and spin polarization in the proposed switching process. b Magnetization dependence on the electric polarization during the switching process. Fitting is marked by solid line. ΔQ gives the value difference of magnetic order parameter between DFT calculation and the fitting. The inset shows a good linear dependence between the magnetic and electrical order parameters.

Full size image

At stage V, the electric polarization reverses its direction to point upward and increases in magnitude, while the magnetic moment of Hf1 (Hf3) starts to decrease (increase) in magnitude due to the similar reason discussed in Fig. 2c and a net magnetization emerges and reverses its direction to the right. At stage VI, with the further decrease of the magnetic moment of Hf1, the net magnetization changes its direction again and points to left. At stage VII, the system goes into an equivalent state in which the electric polarization and net magnetization arrive to the same values but with their directions reversed with respect to their initial states at stage I. Such direct correlation between ferroelectricity and ferrimagnetism suggests an inherent magnetoelectric coupling in the system.

To get an insight of the magnetoelectric coupling in this system, a phenomenological potential is given below based on the behaviors of two coupled order parameters shown in Figs. 4b and S7. As implied above and by Supplementary Fig. 7, it is clear that the net spin magnetization/polarization is stabilized by the electric polarization and thus considered as a secondary order parameter. To simplify the form of the potential, only two order parameters are included. One order parameter \({Q}_{{\varGamma }^{-}}\) denotes the electric polarization and the other one \({Q}_{{mG}}\) characterizes the net spin polarization (i.e., proportional to MHf1 + MHf2 + MHf3 + MHf4). Here, \({Q}_{{\varGamma }^{-}}({Q}_{{mG}})=\pm 1\) corresponds to the electric polarization (spin polarization) of the ground states and different signs mean the reversal of their direction. Neglecting high-order coupling terms and considering the fact that the electric (spin) polarization here is primary (secondary), the potential is given by

$$\begin{array}{ll}{H}_{L}=&{H}_{0}+{{C}_{20}Q}_{{\varGamma }^{-}}^{2}+{{C}_{40}Q}_{{\varGamma }^{-}}^{4}+{C}_{02}{Q}_{{mG}}^{2}\\&+\,{C}_{11}{Q}_{{\varGamma }^{-}}{Q}_{{mG}}+{C}_{21}{Q}_{{\varGamma }^{-}}^{2}{Q}_{{mG}}+{C}_{31}{Q}_{{\varGamma }^{-}}^{3}{Q}_{{mG}}+{C}_{12}{{Q}_{{\varGamma }^{-}}Q}_{{mG}}^{2}+{C}_{22}{Q}_{{\varGamma }^{-}}^{2} {Q}_{{mG}}^{2}\end{array}$$
(1)

where \({{H}}_{0}\) is the potential energy of the fictitious high-symmetric phase (where \({Q}_{{\varGamma }^{-}}\) and \({Q}_{{mG}}\) is zero) and \(C\)’s are coefficients. After minimizing the energy with respect to \({Q}_{{mG}}\), one can get

$${Q}_{{mG}}=-\frac{{C}_{11}{Q}_{{\varGamma }^{-}}+{C}_{21}{Q}_{{\varGamma }^{-}}^{2}+{C}_{31}{Q}_{{\varGamma }^{-}}^{3}}{2{C}_{02}+2{C}_{12}{Q}_{{\varGamma }^{-}}+{2C}_{22}{Q}_{{\varGamma }^{-}}^{2}}$$
(2)

For the system studied here, \({C}_{40}\), \({C}_{02}\) and \({C}_{22}\) are positive, while \({C}_{31}\) and \({C}_{20}\) are negative in accord with Figs. 1d, 4a. Generally, as it can be inferred from Eq. (2), the spin polarization (proportional to \({Q}_{{mG}}\)) decreases in an almost linear fashion with the electric polarization (proportional to \({Q}_{{\varGamma }^{-}}\)) when \({Q}_{{\varGamma }^{-}}\) is large in magnitude, while it evolves as the summation of linear and cubic terms of \({Q}_{{\varGamma }^{-}}\) when \({Q}_{{\varGamma }^{-}}\) is small (Fig. 4b). Quantitatively, as clearly shown in Fig. 4b, Eq. (2) is in good agreement with the data points (black square) by the fitting denoted by red solid line (the fitted parameters are given in Supplementary Table 2). Such behavior is rarely reported in 2D materials. The presence of such linear and cubic coupling terms in spin polarization without the inclusion of spin-orbit coupling indicates a strong magnetoelectric coupling, which is in great contrast to that found in conventional multiferroic materials such as BiFeO3, TbMnO3 and in known 2D multiferroics20,22,23,32,34,44,45. Note that the deviation from the fitting line at the end points is likely due to the instability of Stoner magnetism (see Supplementary Fig. 8), which is another manifestation of strong magnetoelectric coupling. Considering the metallic nature of the material, the minor contribution of pure carrier-mediated magnetoelectric effect cannot be ruled out46. Additionally, the inset of Fig. 4b shows the dependence between magnetic and electrical polarization with the ground-state crystal structure. A good linear relationship is clearly observed within the given window, which implies a giant tunability of magnetization through electric polarization or vice versa. Such properties can be practically useful for functional devices.

Proposed functional devices

In Fig. 5, schematics of two potential applications with bilayer Hf2S are shown, taking advantage of its strong magnetoelectric effects and metallic properties. Figure 5a presents a possible magnetic field sensor. The presence of an external perturbative magnetic field can induce a change in the magnetic order that is strongly coupled to the electric polarization (see the inset of Fig. 4b). Then the change in electric polarization can be detected by charge flow. Figure 5b gives an ultrathin device based on magnetoelectric magnetoresistance. An applied electric field can efficiently tune the magnetic moments of Hf2 on the top layer by manipulating the Stoner instability (see Fig. 2c and its discussion). Such tunability can be continuous. Considering the spin-dependent scattering within the layers47, electrical resistivity varies with the change of the magnetic ordering, which may be further utilized in a similar way as conventional resistors.

Fig. 5: Schematics for possible potential applications.
figure 5

a Magnetic field sensor. b Magnetoelectric magnetoresistance.

Full size image

We also note that it may be challenging to separate bulk electrides into individual layers due to the strong interlayer electrostatic interaction. Given the fact that few-layered electrides in ultrahigh vacuum have been obtained experimentally48,49, we also consider the polar stackings of four, six and eight layers of Hf2S, which exhibits less prominent but qualitatively similar behaviors as those in the bilayer case (see Supplementary Fig. 9 and Supplementary Table 3).

In summary, an unexpectedly large polarization contributed by the delocalized interlayer anionic electrons in bilayer Hf2S electride is identified. Meanwhile, a coexistence of ferroelectricity and metallicity, which are supposed to be mutually exclusive in a single-phase material, is predicted for this system. Besides that, the polar stacking bilayer Hf2S also exhibits a net spin polarization resulting from the uncompensated magnetic moments of Hf atoms located at the two opposite surfaces. Furthermore, the spin magnetization is intrinsically generated and coupled to the electric polarization through a complex dependence corroborated by a phenomenological model. The intriguing combination of electric polarization, metallicity and strong magnetoelectric effects of a novel origin in 2D materials is extremely rare in nature. Similar phenomena are also observed in other systems such as the polar stacking bilayer Hf2Te and Zr2S from our further calculations. Our findings therefore extend the emerging realm of 2D multiferroics and provide a new perspective for searching for strong magnetoelectric effects in ultrathin materials.

Methods

Density functional theory calculations

The calculations were performed within the density functional theory (DFT) by using the projector augmented wave (PAW) method, as implemented in the VASP code50,51,52,53. The revised generalized gradient approximation (GGA) with Perdew-Burke-Ernzerhof (PBE)54 was used to treat the exchange-correlation interactions and tested to give good bulk lattice constants43,55. The energy cutoff for the plane-wave basis set is 400 eV (convergence tests have been shown in Supplementary Fig. 10). The \(k\)-point sampling was performed with a Monkhorst–Pack scheme of a \(27\times 27\times 1\) grid with the Γ point included. The convergence criteria for the Hellmann-Feynman force and total energy were set at 0.005 eV/Å and 10−7 eV, respectively, for structural relaxations. The polarization is calculated by a direct integration of electronic charge density multiplied by the position vector over space within a proper slab shape. Phonon dispersion is calculated with the Phonopy package56,57. A vacuum space of about 18 Å was used to minimize false interactions between periodic images and nonlocal van der Waals interaction between layers is captured by the optB86b-vdW functional58. The correlation effects of 5d electrons in Hf were examined and found not to alter our results qualitatively, as seen from Supplementary Table 1. We thus ignore the Hubbard U correction in our calculations, as done in previous reports43,59. The magnetic ground state is confirmed by Monte Carlo simulation60.