Realization of 2D Multiferroic with Strong Magnetoelectric Coupling by Intercalation: A First-principles High-throughput Prediction

The discovery of novel two-dimensional (2D) multiferroic materials is attractive due to their potential for the realization of information storage and logic devices. Although many approaches have been explored to simultaneously introduce ferromagnetic (FM) and ferroelectric (FE) orders into a 2D material, the resulting systems are often plagued by weak magnetoelectric (ME) coupling or limited room-temperature stability. Here, we present a superlattice strategy to construct non-centrosymmetric AM 2 X 4 multiferroic monolayers, i.e., intercalating transition metal ions (A) into the tetragonal-like vacancies of transition metal dichalcogenide bilayers (MX 2 ). Starting from 960 intercalated AM 2 X 4 compounds, our high-throughput calculations have identi�ed 21 multiferroics with robust magnetic order, large FE polarization, low transition barrier, high FE/ FM transition temperature, and strong ME coupling. According to the origin of magnetism, we have classi�ed them into twelve type-a, seven type-b, and two type-c multiferroics, which also exhibit different ME coupling behavior. During the switching of polarization, the reversal of skyrmions chirality, the transition of magnetic ground state from FM to antiferromagnetic, and the changes in spin polarized electron spatial distribution were observed in type-a, type-b, and type-c 2D multiferroic materials, respectively. These results substantially expand the family of 2D ferroic materials and pave an avenue for designing and implementing nonvolatile logic and memory devices.


Introduction
6][7] They are a special class of materials that simultaneously exhibit two or more primary ferroic orders, such as ferroelectricity, ferromagnetism, ferroelasticity, and ferrovalley, and there exists coupling between them, which can be mutually controlled.In particular, the coupling between ferroelectric (FE) and ferromagnetic (FM) orders leads to the magnetoelectric (ME) coupling effect, which offers the unique opportunity to control the spin by applying an electric eld.10] Depending on the microscopic FE mechanism, multiferroic materials can be divided into two categories, i.e., type-I and type-II.In type-I multiferroics, the ferroelectricity and magnetism originate from different atoms.The coupling between electric dipoles and magnetic spins mediated by spin-lattice interaction usually leading to a weak ME coupling, but often exhibits considerable spontaneous polarization and high transition temperature.To date, several type-I multiferroics such as CrN, 11 M I M II P 2 X 6 (M I , M II = metal elements, X = O/S/Se/Te), 12 and ReWCl 6 13 have been predicted theoretically, while multiferroic has been con rmed experimentally in 2D CuCrP 2 S 6 . 14In type-II multiferroics, FE order generally arises from special magnetic order, and the direct correlation between polarization and magnetization promises strong ME coupling.However, special magnetic order, such as spiral or complex collinear magnetic order, could lead to spin and charge frustration, which signi cantly reduce the intensity of FE polarization and transition temperature.For example, the 120° non-collinear Y-type antiferromagnetic (AFM) order of Hf 2 VC 2 F 2 monolayer can induce polarization perpendicular to the spin spiral plane, approximately 2.9 × 10 − 7 µC/m, much smaller than most 2D ferroelectrics. 15There was also experimental evidence that NiI 2 monolayer exhibited proper screw state with a given handedness coupled to the charge degrees of freedom, producing a chirality-controlled electric polarization with a transition temperature of 21 K. 16 Since both type-I and type-II multiferroics have advantages and shortcomings, continuous search for single-phase 2D multiferroic materials with excellent FE-FM multiferroic properties, i.e., strong ME coupling, high Curie temperature, and FE polarization intensity, is still a challenge for this fast-growing eld.
Intercalation is a chemical process in which foreign species insert into crystal gaps to eventually form hybrid intercalation compounds.In 2D layered materials with van der Waals (vdW) gap, intercalation reaction has been proven as a powerful approach for synthesizing single-phase crystals with the coexistence of multiple order parameters. 17A family of 2D MA 2 Z 4 materials was obtained by intercalating a MoS 2 -type MZ 2 monolayer into InSe-type A 2 Z 2 bilayers.Among them, 2D MnBi 2 Te 4 is an AFM axion/Chern insulator with quantized anomalous Hall conductivity, 18-20 while PdBi 2 Te 4 is a superconducting topological metal. 21Nontrivial topological properties, FM semiconductors, Ising superconductivity, and robust electron valleys are also found in SrGa 2 Te 4 , NbSi 2 N 4 , TiSi 2 N 4 , MoSi 2 P 4 , MoSi 2 As 4 , WSi 2 P 4 , and WSi 2 As 4 . 18,22In addition, the Kagome lattice formed by Ta 8 Se 12 (66.7%Taintercalated) stabilizes the charge-density wave and exhibits FM metallicity, 23 while Cr 1/3 TaS 2 (Crintercalated 2H-TaS 2 ) is a chiral crystal, which has a helical spin texture and continuously transforms into a chiral soliton lattice, eventually reaches a forced FM state under an external magnetic eld. 24In addition, 2D AgCrS 2 composed of Ag intercalated CrS 2 bilayers has been experimentally demonstrated to be stable and exhibit superionic and FM behavior at room temperature. 25,26Brie y speaking, intercalation breaks the inversion symmetry of the system, changes the bonding order from vdW to ionic or covalent type, and modulates the band structures of host lattice and intercalated species.
Inspired by the exotic properties triggered by unique intercalated architectures, herein we propose a general approach to systematically design 2D multiferroic materials by non-centrosymmetric intercalation of 3d/4d transition metal ions (A) into the 2H/1T phase of TMD bilayers (MX    27 Inspired by these results, we have constructed a series of hybrid intercalation compounds, in which the 2H or 1T phase MX 2 bilayers with AB stacking is the host lattice, and A atoms is the guest species.When the intercalated A atoms uniformly occupy the tetragonal-like vacancies in the interlayer of bilayer MX 2 , an A atom coordinates with one X atom from the top layer and three X atoms from the bottom layer.This leads to an unequal distance between the A atom and two MX 2 layers, thereby breaking the central inversion symmetry of the system and generating a spontaneous polarization.The tetrahedral crystal eld can be switched by a 180° rotation when an A ion bonds with three X atoms in the top layer and one X atom in the bottom layer, leading to polarization reversal.The selection of MX 2 layer and A atoms is summarized in Supplementary Table 1, and a total of 960 non-centrosymmetric intercalation compounds were obtained with all combinations of MX 2 layer and A atom.The representative crystal structures of 2H phase compounds AM 2 X 4 (H-AM 2 X 4 ) and 1T phase compounds AM 2 X 4 (T-AM 2 X 4 ), as well as the corresponding FE switching mechanism are shown in Fig. 1.

Results and discussion
Figure 2 schematically illustrates the procedure of high-throughput DFT calculations on the structures and ferroic properties of 2D 2H/1T-AM 2 X 4 systems.Firstly, we performed a full structural optimization to con rm that the intercalated A ions in all 960 AM 2 X 4 compounds deviated from the center position and caused spontaneous polarization.Secondly, we evaluated the dynamic and thermodynamic stabilities of these AM 2 X 4 compounds by computing their phonon dispersion and formation energy.As can be seen in Supplementary Fig. 1, there is no negative frequency throughout the entire Brillouin zone of 104 2D AM 2 X 4 compounds.Among them, the formation energies of 100 systems are negative (Supplementary mainly contributed by the S atoms and the CBMs stem mainly from the intercalated Co atoms (Supplementary Fig. 6).Except for the aforementioned three semiconductors, all other intercalation compounds exhibit metallic or half-metallic behavior, as the 3d/4d orbitals of the intercalated transition metal ions are usually partially occupied and have strong itinerant characteristics.Based on this nding, we will discuss the intrinsic FE behavior of semiconductors and metals/half-metals separately in the following content.
[33] For metals or half-metals, there only exists P out , since P in is annihilated by the presence of free electrons.
Using metallic T-PdZr 2 Se 4 as an example, we discuss the mechanism of coexistence of ferroelectricity and metallicity.The metallic behavior of T-PdZr 2 Se 4 was con rmed by calculating the electronic band structure and projected density of states (PDOS), as displayed in Figs.3a-b.The conduction electrons calculated by Eq. ( 3) are mainly distributed in the top and bottom layers of the system (Fig. 3c), which is consistent with the PDOS results.However, the polarization electrons obtained by Eq. ( 5) are mainly distributed around the intercalated Pd atoms and exhibit an oscillatory behavior (Fig. 3c), which is related to the non-uniform distribution of charge density caused by the bonding difference between the Pd atoms and their top and bottom bonding Se atoms.The difference in the distribution between polarization electrons and conduction electrons prevents the annihilation of P out , where the conduction electrons distributed around the Zr atoms on the lower layer cannot fully screen the vertical polarization originating from the intercalating Pd atom.Therefore, ferroelectricity and metallic behavior can coexist in T-PdZr 2 Se 4 with calculated P out of about 3.10 pC/m.This mechanism is also responsible for the P out of other seventeen 2D non-magnetic FE compounds with metallic or half-metallic behavior, as summarizes in Table 1, with values ranging from 0.43 to 9.61 pC/m, lager than the experimentally observed polarization of the metallic WTe 2 bilayer (0.42 pC/m).For 40 FE materials after screening, we further examined their magnetic ground states by comparing the total energy of FM with various types of AFM and FiM con gurations.Among them, 21 systems exhibit both ferroelectricity and magnetism, including ten FM, nine AFM, and two FiM members (Table 2).It is worth noting that the presence of magnetism does not eliminate out-of-plane FE behavior in these intercalation compounds, even for the metallic systems.For example, in T-CoTi 2 Te 4 , the contributions from top and bottom TiTe 2 layers to the conduction electrons are similar in the spin-up channel, while the conduction electrons in the spin-down channel originate from the top TiTe 2 layer (Figs.3e-f).In contrast, the polarization charge density mainly accumulates around the intercalated Co atoms, as displayed in Fig. 3f.The different origins of the two types of charges ensure the possibility of P out in metallic 2D FE systems, no matter whether the system exhibits magnetism.
Table 2 The band gap (eV) at GGA + U level, ground state (GS), out-of-plane polarization (P out ) per unit cell in pC/m, FE transition barrier (E B ) in meV/u.c., magnetic atom (M-atom), average magnetic moment M (µ B / atom), easy magnetization direction EMA, magnetic anisotropy energy MAE (meV/magnetic atom) and Curie temperature T C (K) or Néel temperature T N (K) simulated by Monte Carlo for 21 multiferroic materials, where Ti top and Ti bot represent the Ti atoms in the top and bottom TiX 2 layers of T-CoTi 2 X (X = Se/Te), respectively.[100] −2.05 −

Material
In the following, we focus on the 21 intercalation materials with coexistence of ferroelectricity and magnetism, and perform a detailed analysis of their magnetism and ME coupling effects.For simplicity, we classify the 21 systems into three types according to the origin of magnetism: (i) type-a, in which the magnetism of the system originates from the magnetic atoms M in the top and bottom The Heisenberg model was adopted to describe the magnetic interaction: where S i and S j are the magnetic moments at i and j sites, respectively, J ij is the magnetic exchange coupling parameter between i and j sites, A z is the single-ion anisotropic energy, and the sum over i and j passes through all the magnetic ions.Considering the possible magnetic exchange interaction between the top and bottom MX 2 layers, exchange interactions up to fth nearest neighbors are involved by comparing the energies of FM and different AFM con gurations (Supplementary Fig. 7a and Supplementary Table 4).Here , and represent the intralayer coupling, while and are the interlayer interactions, as summarized in Supplementary Table 5. Taking T-CdCr 2 Te 4 as a representative, and contribute most signi cantly to FM state, which can be explained by the Goodenough-Kanamori-Anderson (GKA) rule, 35 while the effect of is too weak to affect the intralayer FM order.As for the interlayer magnetic interactions, both and favor FM order.
Therefore, the Cr-Cr pairs in the intralayer and interlayer exhibit FM interactions.The magnetic anisotropy energy (MAE) is also an important parameter of a magnetic system that maintains long-range magnetic order against thermal uctuation.It is de ned as the energy difference between the magnetic con gurations with in-plane and out-of-plane magnetization: MAE = E in−plane − E out−plane .Hence, a positive   Moreover, the breaking of time-reversal symmetry ensures the presence of Heisenberg exchange interactions, while the strong spin-orbit coupling (SOC) induced by heavy Te element and the intrinsic inversion symmetry breaking generated by Cd atoms lead to a signi cant DMI.This naturally leads to the formation of magnetic skyrmions.Therefore, we investigated the DMI of two equivalent FE phases, and explored the possible topological magnetic texture and ME coupling behavior in T-CdCr 2 Te 4 .After considering DMI, the spin Hamiltonian can be rewritten as follows: where the magnetic coupling parameter can be extracted by the four-state energy mapping method, 36 and is the DMI vector.The DMI components between the nearest Cr-Cr pairs of the top and bottom Cr layers have opposite directions (Figs.4a and 4h).Since is a directional vector, we calculated and along x, y, and xy directions for the rst nearest neighbor Cr-Cr pairs of the top and bottom layers, respectively, as shown in Table 3 and Fig. 4. The D value of the top Cr layer in the FE1 state is about three times larger than that of the bottom layer, while it is reversed in the FE2 state, indicating that DMI is very sensitive to structural asymmetry.Moreover, the D/|J| ratios of the top and bottom Cr layers in the FE1 state are about 9.03% and 20.66%, respectively, while they change to 20.80% and 9.52% in the FE2 state.Such a large DMI and D/|J| ratio may give rise to peculiar topological magnetic textures like magnetic skyrmions.At the same time, D/|J| values of the top and bottom Cr layers in the FE1 and FE2 states changes with switching of polarization.Therefore, one expects that FE polarization can regulate the topological magnetic behavior in different Cr layers.4, showing a wealth of topological magnetic excitations.As shown in Fig. 4b, the top Cr layer in FE1 state are exhibit a labyrinth-like spin texture (stripe domains) at low external elds (0 ~ 1.0 T), and the edge of the stripe domains wall becomes blurred with the increase of out-of-plane magnetic eld strength.When the applied magnetic eld reaches 2.4 T, the domain structure is broken and anti-skyrmions appear (Fig. 4c), which can be maintained under a magnetic eld of 2.4 ~ 3.17 T.Moreover, the diameter of anti-skyrmions decreases with the increase of magnetic eld, i.e., from 8.9 nm at 2.4 T to 3.8 nm at 3.17 T, which is much smaller than the reported diameters of 2D Cr(I, Cl) 3 (10.5 nm at 0.8 T), 37 Cr 2 Ge 2 Te 6 (77 nm at 0.2 T), 38 and the typical size of skyrmions in thin lms (10 to 100 nm). 39Compared to domain walls in racetrack memory, the smaller size of skyrmions allows for faster data ow or higher data throughput with the same current density, resulting in higher information density.Therefore, T-CdCr 2 Te 4 stands out as an excellent candidate for low-power consumption spintronic devices.When the magnetic eld exceeds 3.5 T, the anti-skyrmions disappear and the FM state becomes favorable (Fig. 4f).In other words, the spin structure of 2D T-CdCr 2 Te 4 undergoes a phase transition with increasing magnetic eld: stripe domain → skyrmion lattice → FM phase.In contrast to the top Cr layer, the bottom Cr layer does not show topological magnetic phenomena under any magnetic eld, which can be attributed to the difference in DMI between the top and bottom Cr layers.In the FE2 state, the top Cr layer behaves the FM phase, while the bottom layer exhibits similar phase transition process as the top layer in FE1 state, and the anti-skyrmions disappear when the magnetic eld reaches 3.7 T, as illustrated in Figs.4i-n.This indicates that the magnetic behavior of the top and bottom layers is reversed by polarization switching.Remarkably, MX 2 layer shows skyrmions with opposite chirality under different polarization directions (up or down) (Fig. 4), indicating that the FE polarization reversal can not only reversibly generate and annihilate magnetic skyrmions, but also achieve chirality reversal of skyrmions, demonstrating interesting ME coupling phenomena.So far, there have been many studies on the control of magnetic skyrmions through FE polarization from both experimental and theoretical aspects.1][42][43] To the best of our knowledge, there was a theoretical proposal of manipulating the magnetic skyrmions in an intercalation compound CoMo 2 S 4 , 27 similar to the present study.However, the chirality of skyrmions was tailored by FE polarization reversal of the intermediate Co layer, which is still different from our strategy in which the modulation of topological magnetic textures occurs in the top and bottom CrTe 2 layers.In other words, this work provides a new approach to polarization-controlled topological magnetism.
We further investigated the effects of temperature and magnetic eld on the topological spin textures of T-CdCr 2 Te 4 monolayer.The magnetic behavior of T-CdCr 2 Te 4 with temperature and magnetic eld evolution is displayed in Fig. 5.At low temperatures (0 K ≤ T < 60 K), the system maintains a stripe domain phase within a small external magnetic eld and transforms to the skyrmion lattice as magnetic eld increases.It should be noted that the critical strengths of magnetic eld required for the formation and disappearance of the skyrmions decrease with increasing temperature.For example, the antiskyrmions appear at B = 2.1 T and persists until B = 3.17 T at zero temperature, but at T = 25 K, they appear in an external magnetic eld of 1.85 T and transform into FM phase when B > 2.1 T. In contrast, at higher temperatures (T ≥ 60 K), the system exhibits a disordered phase at low magnetic elds and transits into FM phase when magnetic eld is enhanced, indicating that the skyrmion lattice phase can be maintained at around 60 K. Therefore, topological phase transition of 2D T-CdCr 2 Te 4 can be achieved by modulating the external magnetic eld and temperature, thereby generating and annihilating skyrmions.Furthermore, the top Cr layer in FE2 state always maintains an FM state, while the magnetic behavior of the bottom Cr layer is similar to the top Cr layer in FE1 state in terms of temperature and magnetic eld effects, with the critical temperature of magnetic skyrmions of around 60 K (Supplementary Fig. 10).
Among the considered systems, there are seven type-b multiferroic materials, including four T-AM 2 X 4 (CoSc 2 S 4 , CoY 2 S 4 , CoZr 2 S 4 , CoZr 2 Se 4 ) and three H-AM 2 X 4 (CrMo 2 Se 4 , MnMo 2 S 4 , MnMo 2 Se 4 ).In contrast to the type-a multiferroics, the intercalated species in the type-b multiferroics are transition metal atoms (Cr, Co, and Mn) with partially lled d orbitals, which are mainly responsible for the magnetism.The pristine MX 2 (M = Sc, Y, Zr, Mo; X = S, Se) bilayers are non-magnetic, and the intercalation of A (A = Cr, Co, and Mn) atoms causes a transition from non-magnetic to magnetic state in the system.
The preferred magnetic order was determined by comparing the energy of different magnetic con gurations (Supplementary Table 4), and three parameters for magnetic exchange interactions (J 1 , J 2 , and J 3 ) are su cient to describe the magnetic interactions as the magnetism is provided by the intercalated A atom within the same layer.Take T-CoZr 2 S 4 for example, its magnetic ground state is FM with the rst nearest neighbor magnetic coupling coe cient J 1 > 0, which can be explained by the GKA rules.Moreover, J 1 of T-CoZr 2 S 4 is approximately 20 times larger than J 2 and J 3 , which dominates the magnetic interaction between the Co atoms.Meanwhile, the magnetization of T-CoZr 2 S 4 prefers the outof-plane direction with a MAE of 0.61 meV/Co, and the Curie temperature from MC simulation is about 70 K, as illustrated in Supplementary Fig.The ME coupling phenomenon of T-CoZr 2 S 4 was explored by examining the evolution of magnetic behavior during FE polarization reversal process.As shown in Fig. 6, the magnetic ground state and the direction of easy magnetization axis change with the FE switching.In the initial phase of polarization reversal, the energies of FM are smaller than that of AFM

Electron density distribution
The conduction electron density is de ned as the charge density within 0.05eV near the Fermi level, and the conduction electron density distribution along z-direction was obtained by integrate within the xy plane: 3 The existence of center inversion symmetry in the paraelectric (PE) state ensures zero polarization, so the net polarization electron density was de ned as the difference in charge density between the FE phase and the PE phase: 4 where and are the total electron density of the FE phase and the PE phase, respectively.Similar to the conduction electron density, we also integrate the polarization charges within the xy plane to obtain the distribution of the polarization electron density along the z-direction: 5

Computational details
All calculations were performed using the density functional theory (DFT), as implemented in the Vienna ab initio simulation package (VASP). 44The ion-electron interactions were described using the projectoraugmented plane wave (PAW) approach, 45 and the electron exchange-correlation interactions were described within the generalized gradient approximation parameterized by Perdew, Burke and Ernzerhof (PBE). 46The planewave basis was expanded to 500 eV, with a convergence threshold of 10 − 7 eV in Monkhorst-Pack k-point grid, 47 with uniform spacing of 2π*0.02Å −1 .A vacuum space of more than 15 Å in z direction was added to avoid interactions between neighboring layers.The dynamic stability of all structures was assessed by their phonon dispersions, which were computed with the nite displacement method implemented in the Phonopy code. 48The thermodynamic stability was evaluated by the formation energy, which is de ned as E f = (E AM2X4 -E A -2 × E M -4× E X )/n, where E AM2X4 is the total energy of AM 2 X 4 monolayers, E A , E M , and E X are the energies of the individual A, M, and X atoms in their most stable bulk phase, respectively, and n = 7 is the total number of atoms per AM 2 X 4 formula.The FE energy transition barrier was evaluated using the climbing-image-nudged elastic band (CI-NEB) method. 49The out-of-plane polarization was calculated considering the dipole correction 50 and the inplane polarization was estimated using the Berry phase method. 51To account for the strong correlation effects involving localized d electrons, the effective Hubbard parameters U eff = U − J were applied to all 3d and 4d transition metal elements, where U is the Coulomb repulsion and J is the exchange parameter.52-54 The U eff values used in this work are summarized in Supplementary Table 1.The Curie or Néel temperature was evaluated using the Monte Carlo (MC) method 55 with a 80×80 lattice.The FE switching temperatures were estimated by the ab initio molecular dynamics (AIMD) simulations within the NVT ensemble, which were performed on a 4×4×1 supercell for at least 10 ps with a time step of 1 fs. 56

Micromagnetic simulations
As for the topological magnetism, we performed MC simulations for CdCr 2 Te 4 monolayer under different temperatures and magnetic elds to simulate the evolution of the spin textures using the Metropolis algorithm based Spirit software. 57The MC simulation considers periodic boundary conditions and operates in an 80×80 supercell containing 115200 spin sites with a minimum MC step size of 6×10 5 , and the initial magnetic moments were set to random con gurations.

34
MX 2 layers; (ii) type-b, in which the magnetism originates from the intercalated A atoms; (iii) type-c, in which both the M atoms in the MX 2 layers and the intercalated A atoms provide magnetism simultaneously.Type-a multiferroic materials, including T-PdCr 2 Te 4 , T-CuMn 2 Se 4 , T-AgMn 2 S 4 , T-AgMn 2 Se 4 , T-AgTc 2 S 4 , T-AgTc 2 Se 4 , T-AgCr 2 Te 4 , T-ZnCr 2 Se 4 , T-CdMn 2 Se 4 , T-CdCr 2 S 4 , T-CdCr 2 Se 4 , and T-CdCr 2 Te 4 , are obtained by inserting Pd/Cu/Ag/Zn/Cd atoms into T-CrX 2 , T-MnX 2 , or T-TcX 2 bilayers.The intercalated A (A = Pd, Cu,Zn, Ag, and Cd) atoms with almost completely occupied 3d/4d orbitals do not contribute to magnetism, while the 3d orbitals of M (M = Cr, Mn and Tc) atoms are partially lled and have more localized electrons and larger magnetic moments.Therefore, the magnetism comes from the M atoms in the MX 2 bilayers, and the special bonding of the intercalated A atom results in negligible differences in the magnetic moments of the M atoms in the top and bottom layers (~ 0.1 µ B ).

(
negative) MAE represents an out-of-plane (in-plane) easy magnetization axis.The calculated MAE of T-CdCr 2 Te 4 is − 0.34 meV/Cr, indicating an in-plane magnetization orientation.Based on the magnetic coupling parameters and MAE, we performed MC simulations with Heisenberg model to estimate the Curie temperature.As displayed in Supplementary Fig. 8a, T-CdCr 2 Te 4 shows intrinsic FM order with a Curie temperature of 260 K, close to room temperature.AIMD simulations for T-CdCr 2 Te 4 at different temperatures was performed to examine the FE stability by analyzing displacement of the intercalated Cd atoms relative to the centrosymmetric position along out-of-plane direction.As shown in Supplementary Fig. 8b, the average displacement of the Cd atoms ( Cd ) in z-direction follows a Gaussian distribution.At T = 250 K, Cd = 0.376 Å represents a macroscopic FE phase.When the temperature reaches T = 300 and 350 K, Cd rapidly decreases to 0.020 and 0.013 Å, respectively, and the peak position gradually moves toward zero and approaches the non-polar phase.These results indicate that the FM and FE transition temperatures of 2D T-CdCr 2 Te 4 are comparable to room temperature, exhibiting good multiferroic stability.In addition, T-AgMn 2 S 4 , T-AgMn 2 Se 4 , T-ZnCr 2 Se 4 , T-CdCr 2 S 4 , T-CdCr 2 Se 4 , T-CdMn 2 Se 4 , and T-CuMn 2 Se 4 also belong to the FM systems with Curie temperatures above room temperature, and T-AgMn 2 Se 4 has the highest Curie temperature up to 525 K ( 9j.Moreover, the magnetic ground state of T-CoZr 2 Se 4 is also FM, while T-CoSc 2 S 4 , T-CoY 2 S 4 , H-CrMo 2 Se 4 , H-MnMo 2 S 4 , and H-MnMo 2 Se 4 favor AFM order, and T-CoY 2 S 4 has the highest Néel temperature of 50 K (Table 2).Except for H-MnMo 2 Se 4 , most structures belonging to type-b multiferroics prefer the out-of-plane direction as the easy magnetization axis, and T-CoY 2 S 4 has a maximum MAE of 0.92 meV/Co.
polarization direction may affect the spatial distribution of spin-polarized electrons.The magnetic ground state of T-CoTi 2 Te 4 maintains FiM order in FE1 and FE2 phases, but the easy magnetization axis along the x-direction in FE1 state transforms to the z-direction after polarization reversal, further con rming the characteristics of ME coupling.To summarize, we have performed a high-throughput rst-principles calculation to lter potential 2D multiferroics induced by displacement of intercalated ion among TMD intercalation compounds.By intercalating transition metal ions into the tetragonal-like vacancies of transition metal dichalcogenide MX 2 bilayers, we obtained 960 candidate 2D structures of AM 2 X 4 , and identi ed40 systems possessing high structural stability, large FE polarization, low FE switching barrier, and high transition temperature.The intrinsic stabilization mechanism of polarization in four FE semiconductors and 36 FE metals or half-metals were carefully discussed and their values of out-plane polarization vary from 0.26 pC/m to 15.22 pC/m, comparable to the typical 2D FE materials.Among 40 FE intercalation compounds, there are 21 materials exhibit ferroelectricity and magnetism simultaneously.According to the origin of magnetism, we divided them into twelve type-a, seven type-b, and two type-c multiferroic with different ME coupling behavior.In type-a multiferroic T-CdCr 2 Te 4 , the magnetism mainly comes from the CrTe 2 bilayer and the conversion of skyrmions between the top and bottom CrTe 2 layers can be achieved by polarization switching.In type-b multiferroic T-CoZr 2 S 4 , the un lled d orbitals of the intercalated Co atoms mainly contribute to the magnetic moment.During the FE switching process, the magnetic ground state of T-CoZr 2 S 4 undergoes a transition from FM to AFM to FM, and the easy magnetization axis undergoes a switch between out-of-plane and in-plane directions.The intercalated Co atoms are able to introduce magnetism into the non-magnetic TiSe 2 and TiTe 2 bilayers, forming type-c multiferroics.In type-c multiferroic T-CoTi 2 Te 4 , both the spatial distribution of the spin-polarized electrons and the MAE change during the polarization reversal.All these results signi cantly expand the family of 2D ferroic materials and open a pathway for the development and implementation of nonvolatile logic and memory devices.
energy and 10 − 3 eV/Å in force, respectively.The Brillouin zone was sampled with the Γ-centered

(
a) and (b) are schematic diagrams of the A atom intercalated into the 1T-and 2H-TMD (MX 2 ) bilayers,where A and A′ represent the atomic positions of the intercalated atoms in FE1 and FM2 states, respectively.In FE1 state, the A atoms form bonds with two X atoms of the bottom MX 2 layer and one X atom of the top MX 2 layer, while in FE2 state, the A′ atom is bonded to two X atoms of the top MX 2 layer and one X atom of the bottom MX 2 layer.The blue arrows describe the relative displacement of the intercalated atoms between FE1 and FE2 states.

Figure 4 The
Figure 4

Figure 5 The
Figure 5

Figure 6
Figure 6 2 ).By ltering a total of 960 possible combinations, we have reproduced the previously reported CoMo 2 S 4 , 27 AgCr 2 S 4 and AgCr 2 Se 4 ,26and further unveiled 40 dynamically and thermodynamically stable 2D AM 2 X 4 ferroelectrics with low FE switching barrier (< 200 meV/f.u.).Among them, different magnetic orders such as FM, AFM and ferrimagnetic (FiM) are identi ed in the 21 FE systems.Depending on the position of the magnetic ions (intercalated A ions or host M atoms), we classify these 21 multiferroics into three types and investigate their different ME coupling behavior.Our intercalation strategy opens a new avenue to develop novel 2D materials with multiferroic properties and explore other fascinating topological phenomena.

Table 1
The band gap (eV) at GGA + U level, ground state (GS), spontaneous out-of-plane polarization (P out ) per unit cell in pC/m, FE transition barrier (E B ) in meV/u.c. of 19 non-magnetic FE intercalation compounds.

Table 2
easy magnetization axis, as summarized in Table 2, T-CdMn 2 Se 4 has the largest MAE of − 0.99 meV/Mn among FM materials, the MAE of T-AgTc 2 Se 4 with AFM order up to − 2.49 meV/Tc, while T-AgTc 2 S 4 with in-plane easy magnetization axis has an MAE of up to 2.09 meV/Tc.
states, meaning that T-CoZr 2 S 4 keeps FM ground state.However, as it approaches the PE phase, the energies of AFM states decrease, and the magnetic ground state becomes AFM state.Afterwards, the energy difference between AFM and FM states increases, and the magnetic ground state returns to FM state again.Interestingly, the change in MAE during polarization reversal process exhibits similar behavior, that is, T-CoZr 2 S 4 favors in-plane MAE with a larger value of − 13.34 meV/Co in PE phase, but adopts an out-of-plane easy magnetization axis in other states of the polarization reversal process.The magnetic ground state and MAE values for the other type-b multiferroics are also summarized in Supplementary Table 5.During the polarization reversal process, H-CoY 2 S 4 and H-MnMo 2 S 4 exhibit similar ME coupling phenomenon as T-CoZr 2 S 4 , while H-CoSc 2 S 4 only changes its magnetic ground state, and H-CoZr 2 Se 4 and H-MnMo 2 Se 4 only change the magnetic anisotropy (Supplementary Fig.11).In type-c multiferroics, pristine TiX 2 (X = Se/Te) bilayers are non-magnetic semiconductor, and the intercalation of magnetic Co atoms produces two stable multiferroic materials with lower barriers, namely, T-CoTi 2 Se 4 and T-CoTi 2 Te 4 .In addition to the intercalated Co atoms with intrinsic magnetic moment, there are also induced moments on the Ti atoms, and the special coordination between the intercalated Co atoms and the Te atoms results in distinct magnetic moments between the Ti atoms in top (Ti top ) and bottom (Ti bot ) layer.Therefore, one has to consider possible FiM con gurations.The magnetic con gurations and corresponding energies of FM, AFM and FiM states for type-c multiferroics are presented in Supplementary Fig.7cand Supplementary Table6.Both T-CoTi 2 Se 4 and T-CoTi 2 Te 4 exhibit FiM ground states with total net magnetic moment of 0.24 µ B /f.u. and 0.21 µ B /f.u., respectively.Using T-CoTi 2 Te 4 as a representative, we have investigated the possible ME coupling phenomena in typec multiferroic material.The electronic band structures of FE1 and FE2 states exhibit metallic properties (Figs.7a-b), but the conduction electrons at the Fermi level of FE1 state are mainly contributed by spindown electrons, whereas the electrons near the Fermi level of FE2 state are mainly dominated by spin-up channel, as illustrated in Figs.7c-d.Moreover, the magnetic moments of Ti top and Ti bot are 0.62 and 0.82 µ B in FE1 state, whereas they become 0.82 and 0.62 µ B in FE2 state, indicating that the reversal of