# Fileset

[s41467-024-47373-7.pdf](https://mdr.nims.go.jp/filesets/baefb6a4-8c1a-4366-a194-913ffa92e379/download)

## Creator

Qifeng Hu, Yuqiang Huang, Yang Wang, Sujuan Ding, Minjie Zhang, Chenqiang Hua, Linjun Li, Xiangfan Xu, Jinbo Yang, Shengjun Yuan, [Kenji Watanabe](https://orcid.org/0000-0003-3701-8119), [Takashi Taniguchi](https://orcid.org/0000-0002-1467-3105), Yunhao Lu, Chuanhong Jin, Dawei Wang, Yi Zheng

## Rights

[Creative Commons BY Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0/)

## Other metadata

[Ferrielectricity controlled widely-tunable magnetoelectric coupling in van der Waals multiferroics](https://mdr.nims.go.jp/datasets/1be6df20-93ec-4700-9b58-550ac1f97b32)

## Fulltext

Ferrielectricity controlled widely-tunable magnetoelectric coupling in van der Waals multiferroicsArticle https://doi.org/10.1038/s41467-024-47373-7Ferrielectricity controlled widely-tunablemagnetoelectric coupling in van der WaalsmultiferroicsQifeng Hu 1,9, Yuqiang Huang1,9, Yang Wang 2,9, Sujuan Ding3,9,Minjie Zhang1,9, Chenqiang Hua1, Linjun Li 4 , Xiangfan Xu 5, Jinbo Yang 6,Shengjun Yuan 7, Kenji Watanabe 8, Takashi Taniguchi 8, Yunhao Lu 1 ,Chuanhong Jin 3 , Dawei Wang 2 & Yi Zheng 1The discovery of various primary ferroic phases in atomically-thin van derWaals crystals have created a new two-dimensional wonderland for exploringand manipulating exotic quantum phases. It may also bring technical break-throughs in device applications, as evident by prototypical functionalities ofgiant tunneling magnetoresistance, gate-tunable ferromagnetism and non-volatile ferroelectric memory etc. However, two-dimensional multiferroicswith effective magnetoelectric coupling, which ultimately decides the futureofmultiferroic-based information technology, has not been realized yet. Here,we show that an unconventional magnetoelectric coupling mechanism inter-locked with heterogeneous ferrielectric transitions emerges at the two-dimensional limit in van der Waals multiferroic CuCrP2S6 with inherent anti-ferromagnetism and antiferroelectricity. Distinct from the homogeneousantiferroelectric bulk, thin-layer CuCrP2S6 under external electric field makeslayer-dependent heterogeneous ferrielectric transitions, minimizing thedepolarization effect introduced by the rearrangements of Cu+ ions within theferromagnetic van der Waals cages of CrS6 and P2S6 octahedrons. Theresulting ferrielectric phases are characterized by substantially reducedinterlayer magnetic coupling energy of nearly 50% with a moderate electricfield of 0.3 V nm−1, producing widely-tunable magnetoelectric coupling whichcan be further engineered by asymmetrical electrode work functions.The interplay and mutual control of multiple ferroic orders, in parti-cularly magnetism and ferroelectricity, in a single-phase system pro-vides a fascinating platform for both fundamental understanding oncorrelated phenomena and developing the next-generation non-vola-tile storage technologies1. Despite the tremendous success in classi-fying and revealing multiferroic mechanisms in paradigmatic bulksystems, viz. type-I multiferroics with nearly independent magnetismand ferroelectricity2 and type-II multiferroics originated fromDzyaloshinskii-Moriya interaction and other types of magnetic order-ings such as collinear magnetic structures3–5, the quest for a room-temperature multiferroic with efficient coupling between magnetismand ferroelectricity remains as the holy grail of the resurgent field.Recent studies on magnetism and ferroelectricity in two-dimensional(2D) van der Waals (vdW) crystals open unparalleled opportunities formultiferroic research6–10, considering that abundant choices of ferroicatomic layers with single-crystal quality can be readily tuned by elec-trostatic field, strain, interfaces and spin-lattice interactions11–14.Equally important, the Lego combination of different 2D ferroicbuilding blocks may lead to the finding of novel multiferroic phasesbeyond the prevailing spin-driven ferroelectric paradigm1,15,16.Received: 26 July 2023Accepted: 27 March 2024Check for updatesA full list of affiliations appears at the end of the paper. e-mail: lilinjun@zju.edu.cn; luyh@zju.edu.cn; chhjin@zju.edu.cn; phyzhengyi@zju.edu.cnNature Communications |         (2024) 15:3029 11234567890():,;1234567890():,;http://orcid.org/0000-0002-2824-5398http://orcid.org/0000-0002-2824-5398http://orcid.org/0000-0002-2824-5398http://orcid.org/0000-0002-2824-5398http://orcid.org/0000-0002-2824-5398http://orcid.org/0009-0006-2588-7922http://orcid.org/0009-0006-2588-7922http://orcid.org/0009-0006-2588-7922http://orcid.org/0009-0006-2588-7922http://orcid.org/0009-0006-2588-7922http://orcid.org/0000-0002-2734-0414http://orcid.org/0000-0002-2734-0414http://orcid.org/0000-0002-2734-0414http://orcid.org/0000-0002-2734-0414http://orcid.org/0000-0002-2734-0414http://orcid.org/0000-0001-7163-4957http://orcid.org/0000-0001-7163-4957http://orcid.org/0000-0001-7163-4957http://orcid.org/0000-0001-7163-4957http://orcid.org/0000-0001-7163-4957http://orcid.org/0000-0003-3517-9701http://orcid.org/0000-0003-3517-9701http://orcid.org/0000-0003-3517-9701http://orcid.org/0000-0003-3517-9701http://orcid.org/0000-0003-3517-9701http://orcid.org/0000-0001-6208-1502http://orcid.org/0000-0001-6208-1502http://orcid.org/0000-0001-6208-1502http://orcid.org/0000-0001-6208-1502http://orcid.org/0000-0001-6208-1502http://orcid.org/0000-0003-3701-8119http://orcid.org/0000-0003-3701-8119http://orcid.org/0000-0003-3701-8119http://orcid.org/0000-0003-3701-8119http://orcid.org/0000-0003-3701-8119http://orcid.org/0000-0002-1467-3105http://orcid.org/0000-0002-1467-3105http://orcid.org/0000-0002-1467-3105http://orcid.org/0000-0002-1467-3105http://orcid.org/0000-0002-1467-3105http://orcid.org/0000-0001-6825-7206http://orcid.org/0000-0001-6825-7206http://orcid.org/0000-0001-6825-7206http://orcid.org/0000-0001-6825-7206http://orcid.org/0000-0001-6825-7206http://orcid.org/0000-0001-8845-5664http://orcid.org/0000-0001-8845-5664http://orcid.org/0000-0001-8845-5664http://orcid.org/0000-0001-8845-5664http://orcid.org/0000-0001-8845-5664http://orcid.org/0000-0003-3336-8568http://orcid.org/0000-0003-3336-8568http://orcid.org/0000-0003-3336-8568http://orcid.org/0000-0003-3336-8568http://orcid.org/0000-0003-3336-8568http://orcid.org/0000-0002-9487-1151http://orcid.org/0000-0002-9487-1151http://orcid.org/0000-0002-9487-1151http://orcid.org/0000-0002-9487-1151http://orcid.org/0000-0002-9487-1151http://crossmark.crossref.org/dialog/?doi=10.1038/s41467-024-47373-7&domain=pdfhttp://crossmark.crossref.org/dialog/?doi=10.1038/s41467-024-47373-7&domain=pdfhttp://crossmark.crossref.org/dialog/?doi=10.1038/s41467-024-47373-7&domain=pdfhttp://crossmark.crossref.org/dialog/?doi=10.1038/s41467-024-47373-7&domain=pdfmailto:lilinjun@zju.edu.cnmailto:luyh@zju.edu.cnmailto:chhjin@zju.edu.cnmailto:phyzhengyi@zju.edu.cnBeyond the frameworkof the isotropicHeisenbergmodel, variousvdWbulkmagnetswith significantmagnetocrystalline anisotropy havebeen re-explored at the 2D limitwith fruitful results, e.g., 2DCrI3, CrCl3,CrSBr, and CrPS4 with interlayer antiferromagnetic (AFM) couplingand intralayer ferromagnetic (FM) ordering, layer-dependent FM inCr2Ge2Te6, CrBr3 and Fe3GeTe2, and emergent intralayer AFM systemsof MnPS3 and CrOCl13,14,17–27. Gate-tunable magnetic properties havealso been demonstrated by electric controlled magnetic transistors of2D Cr2Ge2Te6 and memristive devices of CrI311,28. As an electric coun-terpart to magnetism, 2D ferroelectricity in atomically thin vdW crys-tals is also not rare, such as SnTe, CuInP2S6, In2Se3, WTe2 and bilayerboron nitrides with divergent ferroelectric mechanisms8,29–34. Intrigu-ingly, unlike bulk materials, 2D vdW systems tend to host multipleferroic orders in one single phase, e.g., ferroelectricity and ferroelas-ticity predicted in monolayer group IV monochalcogenides35,36, andferromagnetism and ferroelasticity in monolayer α-SnO37. Recently,the experimental implementation of a 2Dmagnetoelectricmultiferroicwas reported in NiI29, however, the electrical control of magnetism,and the vice versa, remains elusive in this system.Belonging to the large family of metal thio- and selenophosphatevdW crystals with rich ferroic phases21,29,38,39, CuCrP2S6 (CCPS) is aunique choice of magnetoelectric multiferroics for exploring 2Dmagnetism-ferroelectricity interactions with experimentally con-firmed coexistence of interlayer AFM and intralayer stripe antiferro-electric (AFE) in the bulk crystals40–51. As illustrated in Fig. 1a, CCPScrystals have a monoclinic lattice of the space group Pc (No. 7), inwhich a monolayer can be viewed as hexagonal vdW cages consistingof interconnected CrS6 and P2S6 octahedrons. Below the Néel tem-perature (TN) of ~ 31 K, the magnetic moments of Cr3+ become ferro-magnetically aligned within individual monolayers, which are furtherAFM coupled in perpendicular to the vdW plane (Fig. 1b). On the otherhand, the ferroelectric property of CCPS is rooted in Cu+ ions, whichare randomly distributed within the FM CrS6-P2S6 cages. At 145K, thecrystal becomes AFE, when the adjacent Cu ions form striped AFEalignment driven by a double-well pseudo potential41. In the AFE state,each CCPS monolayer is characterized by a distinctive mirror sym-metry along the b-axis, in perpendicular to the Cu+ AFE stripe chains.Here, we report an unconventional magnetoelectric coupling(MEC) interlocked with electric field tunable ferrielectric (FiE) transi-tions in 2D multiferroic CCPS, unveiling the rich promises of dis-covering highly tunable 2D MEC mechanisms fundamentally differentfrom the bulk counterparts. Using magnetic tunneling junctions(MTJs), we demonstrate the electric field control of interlayermagneticcoupling strength, when Cu+ ions rearrange within the CrS6-P2S6 cagesto form a layer-dependent heterogeneous FiE phase. The existence andcontinuous control of the FiE state are unambiguously proved bytemperature- and polarization-dependent second harmonic genera-tion (SHG) experiments, in which the latter unveils a fingerprinting six-fold to two-fold symmetry transition when thin-layer CCPS is drivenfrom the AFE state into the FiE phase by an external bias voltage (Vb).The formationof the FiE phase enforces substantially lattice distortionsfor the CrS6 octahedrons, which effectively modulate the interlayersuperexchange coupling strength via the interlayerCr-S-S-Cr pathways.This unique FiE-controlled MEC mechanism has a wide tunability, evi-dent by a remarkable reduction in the interlayer magnetic couplingenergy of nearly 50% with a moderate electric field of 0.3 V nm−1.ResultsDue to the insulating nature, we study the layer-dependent 2D mag-netism, ferroelectricity andmagnetoelectric properties of CCPS basedFig. 1 | Magnetoelectric multiferroics and Cu+ ion polarizations in a ferro-magnetic vdW cage. a Top and side views of the lattice structure of a CCPSmonolayer. Interconnected CrS6 and P2S6 octahedrons form hexagonal vdW cagesfor polarizable Cu+ ions. b Schematic of the interlayer AFM order in CCPS, whereonly the Cr atomic layers are shown for clarity. cNCAFM imaging and line profilingof bilayer CCPS flakes on SiO2/Si substrates. d ADF-STEM image of a BL CCPS,showing high crystal quality without noticeable defects. e Top: device structure ofthin-layer CCPS MTJs. Bottom: false-color optical image of a representative BLMTJBL-S44. fTunneling conductanceG vs T curves of representative BL-S44 and 7L-S21.Inset: the corresponding dG/dT, revealing a thickness-independent AFM phasetransition at about 30 K.Article https://doi.org/10.1038/s41467-024-47373-7Nature Communications |         (2024) 15:3029 2onverticalmagnetic tunnel junctions,whichare consistingof few-layergraphene (FLG) electrodes and thin layers of CCPS channels withconsecutive thicknesses from bilayer (BL) to 10L (see Methods andSupplementary Information Note 1). Figure 1c shows the non-contactatomic forcemicroscopy (NCAFM) of a BL CCPS sample, as evident bya 2 nm step height. Atomic-resolution annular dark-field scanningtransmission electron microscopy (ADF-STEM) confirms the highcrystal quality of exfoliated thin-layer CCPS samples with no apparentvacancies (Fig. 1d and Supplementary Note 1). If not specified other-wise, all tunneling results were measured at a base temperatureof T = 1.6K.To investigate the magnetic phase transitions at the 2D limit, wehave carried out T-dependent zero-bias tunneling experiments. Asshown in Fig. 1e, f, using a small excitation current of 10 nA and zeromagnetic field (H), the tunneling conductance of CCPS MTJs with dif-ferent thicknesses consistently show an anomaly at ~ 30K, in excellentagreement with the bulk TN. The results strongly suggest that atom-ically thin CCPS retain the bulk AFM order, as expected for a vdWmagnet with strong magnetic anisotropy. Using H-dependent tunnel-ing spectroscopy, we observe the typical metamagnetic transitionbehavior for thin-layer CCPS MTJs, when interlayer Cr spin momentsbecome symmetrically tilted along the H direction to form the cantedAFM (CAFM) state14,24. As shown in Fig. 2a, the tunneling current (It) ofBL-S44 first increases smoothly as a function of in-plane H∥, mani-festing the spin-canting angle vsH evolution of the CAFMphase. In theAFM ground state, two CCPS MLs with opposite magnetization direc-tion form a spin filter tunneling barrier equivalent for electron trans-missionwith opposite spins, as schemed in Fig. 2b (see SupplementaryNote 2 for detailed information on the spin-filter mechanism). Thisproduces a low It which corresponds to a high resistance state52–54.When μ0H∥ reaches 3.2 T, the atomic spin valve is switched to the lowresistance state due to parallel magnetization between two adjacentMLs. The resulting It becomes significantly enhanced attributed toasymmetric tunneling barriers for opposite spins24,53. For simplicity, wewill refer this magnetization saturation state as the field-aligned fer-romagnetic phase (FM0), and the saturation magnetic field as Hsat.Such an H-driven CAFM to FM0 transition is further confirmed byT-dependent tunneling magnetoresistance (TMR), defined byTMR= [It(H) − It(0)]/It(0) × 100%14,24,53. As summarized in Fig. 2c, theCAFM-FM0 transitions gradually shift to lower values of Hsat due toincreased thermal fluctuation. Above TN, the TMR curves exhibit anegative growth rate, in consistent with a paramagnetic state. It isilluminating to plot the field derivative of TMR (dR/dH∥) into a 2Dcontour image as a function of both T andH∥, which clearly reveals themagnetic phase diagram of CCPS at the 2D limit (Fig. 2d). For CCPSMTJs with odd-layer thicknesses, for example a quintuple-layer MTJ(5L-S5), the corresponding T-dependent TMR and the 2D contour plotare shown in Fig. 2e, f, respectively. The overall phase diagram isagreeing with the BL device, indicating the same CAFM to FM0 phasetransitions. The distinctive odd-layer effect below 2 T is due to a field-induced perfect collinear AFM state which is first reported in few-layerCrCl3 MTJs24. Most critically, for an easy-plane magnetic anisotropydominated vdW AFM, Hsat is directly proportional to the interlayerexchange coupling energy J⊥24. Compared to CrCl3 with the samethickness,Hsat of CCPS is nearly three times higher, which is crucial forthe observation of the FiE-interlocked MEC mechanism.Remarkably, CCPS MTJs show Vb polarity-dependent Hsat, whichreaches a prominent changeof 50%with amoderate electricfield (E) of0.3 V nm−1 for octuple-layer CCPS. As shown in Fig. 3a for MTJ 8L-S20,thenormalizedTMRcurves exhibit anHsat value of 5.9T forVb = − 1.2 V,which drastically changes to 8.8 T for Vb = 1.2 V (see SupplementaryNote 3 for TMR normalization procedures). Distinct from aFig. 2 | Layer-dependent metamagnetism and even-odd effects in thin-layer CCPS. a Tunneling current It vs in-plane magnetic fieldH∥ characteristics of arepresentative bilayer MTJ (BL-S44). BL-S44 exhibits a CAFM saturation field of 3.3T, which is typical for BL devices. Here, green and orange arrows indicate the H∥ramping directions. b Illustration of the multiple spin-filter effect in vdW magnetswith interlayer AFM ordering. c T-dependent TMR vs H∥ of BL-S44, showing thedwindlingof theCAFMsaturationfieldwhenapproachingTN.d2Dcontourplottingof the first derivative of TMR (dR/dH) vs H∥ and T for BL-S44. By taking the deri-vative, the magnetic phase diagram of BL CCPS is clearly revealed. e T-dependentTMR vs H∥ of quintuple-layer MTJ (5L-S12), showing the same dwindling of Hsatwhen approaching TN. f 2D contour of dR/dH vs H∥ and T for 5L-S12. Below 2 T, themagnetic phase diagram of odd-layer CCPS is characterized by a low-field perfectcollinear AFM zone, i.e., the odd-layer effect.Article https://doi.org/10.1038/s41467-024-47373-7Nature Communications |         (2024) 15:3029 3paraelectric dielectric response, the Vb-dependence is conspicuouslynon-linear, which can be divided into two opposite-polarity regimes of0.73 to 1.3 V for highH +sat = 8:8 T and −0.9 to −1.3 V for lowH�sat = 5:9 T,respectively. The dichotomic polarity behavior in Vb and the widetunability of Hsat are best visualized by a 2D contour plot of TMR vsboth Vb and H. As shown in Fig. 3b, the large Hsat window betweenpositive and negative Vb polarities opens a new possibility to directlymanipulate the magnetic state of CCPS MTJs by Vb, i.e., the electricalcontrol ofmagnetism. Indeed, by sweeping Vb andH to follow the pathindicated by red solid lines in Fig. 3b, we can readily enforce a FM0 toCAFM phase transition, and the vice versa, simply by reversing the Vbpolarities. Since Hsat represents the interlayer coupling energy J⊥, wecan directly measure the change of MEC energy introduced by Vb-polarity reversal by ΔEMEC / ðH +sat � H�satÞ. Detailed TMR study revealsthat in between H ±sat, Hsat is continuously tuned by Vb setpoints beforereaching the binary saturation values (Fig. 3c and detailed discussionsin Supplementary Note 4). The large ΔEMEC of ~ 1.08meV per unit-cell(UC) explains the bistable magnetoresistance states TMR± for H ≠0(Fig. 3d). However, within the same AFM ground state, zero-H resis-tance states R ±0 remain binary valued, which can be reproduciblyswitched by reversing the Vb polarities, as shown in Fig. 3e for MTJ10L-S39.Such an unexpected reversal behavior in R ±0 reminds us theimportance of the inherent stripe-AFE phase of CCPS, in analogy withthe resistance switching of graphene-ferroelectric field-effect transis-tors (GFeFETs) under a constant electrostatic doping55,56. However,unlike the hysteretic two-state switching in GFeFETs controlled byferroelectric polarization reversals55,56, the stripe-AFE phase of CCPSrequires an enormous energy to be fully polarized, which can beevaluated quantitatively by the first-principle density-functional the-ory (DFT) calculations. As shown in Fig. 4a, an FE configuration for a BLUC will unrealistically increases the system energy by 337meV per UCwhen compared with the stripe-AFE ground state (see Method for thecalculation details). From the thermodynamic point of view, the FEstate is also not favored due to themaximization of the depolarizationeffect, as depicted by the black hollow arrows. Based on DFT calcula-tions, it is also not feasible to locally flip an AFE stripe by E within asingle vdW ML. For BL CCPS, the lowest metastable state is to flip theanti-parallel stripes of the bottom ML toward the vdW interface, cor-responding to a largely reduced energy requirement of 95meV per UC(2L-mAFE-I). In Fig. 4b, we plot the energy potential diagram of BLCCPS as a function of the displacement (Δd) of the anti-parallel Cu-ionstripes within the bottomML, suggesting a rapid growth in the systemenergy by approaching the Cr atomic plane. Note that even with thepresence of an impractical E of 2 V nm−1, the energy of 2L-mAFE-Iremains substantially higher than the AFE ground state, and the energybarrier in-between is far too high to be crossed by Cu+ ions (Supple-mentary Note 5).Fig. 3 | Ultra-wideVb-tunability ofmagnetoelectric coupling in thin-layerCCPS.a Normalized TMR vs H∥ curves of MTJ 8L-S20 under different Vb setpoints. b 2Dcontour plot of the normalized TMR vs Vb and H∥ in (a). It is noteworthy that foreach H∥ setpoint, the CAFM canting angle θ± for opposite-polarity Vb are divergent(see Supplementary Fig. 5 for CAFM θ± calculations). c ΔEMEC(E) vs Vb, extractedfrom the TMR data of 10L-S39, and the error bars are defined as the ± 10% variationof half-maximum points in dG/dH curve (see Supplementary Figs. 5 and 6).d Dichotomic It vs H∥ characteristics of 8L-S20. Note that It is symmetric on H∥ramping, excluding any ferroelectric polarization reversal. eZero-HR ±0 reversalsbyswitching the Vb polarities. It is apparent that the It −Vb sweeping is completely freeof hysteresis.Article https://doi.org/10.1038/s41467-024-47373-7Nature Communications |         (2024) 15:3029 4Indeed, for all CCPSMTJs, the It −Vb characteristicsmeasuredwithvarious H setpoints (Supplementary Note 6) are thoroughly free ofhysteresis behaviors, in consistent with the DFT calculations thatpredict no metastable phase transition under realistic E values. TheDFT-based analyses of flipping a minimal AFE stripe can be readilyextended to different ML numbers, as shown in the bottom panels ofFig. 4a for trilayer CCPS (quadruple-layer in Supplementary Note 5).With increased layer numbers, different metastable states can be usedas an ML index of the flipped AFE stripe, whose energy cost is closelyrelated to a peculiar layer-dependent depolarization mechanism.Using BL CCPS as an example, by flipping one AFE stripe in the topML,the energy requirement of 2L-mAFE-II is drastically increased by nearlyfour times to 250meV per UC, when compared with 2L-mAFE-I. Thelarge energy difference between 2L-mAFE-I and 2L-mAFE-II is due tothe fact that, for the former, the local flipped AFE domain is effectivelyscreened by the neighboring AFE ML to minimize the depolarizationfield. Such a unique dipole screening mechanism for vdW ferro-electrics, designated as the vdWdepolarization effect hereafter, can bewell reproduced by DFT calculations of layer-dependent charge den-sity distribution in response to AFE stripe displacements (see Supple-mentary Note 5 for details).Although the metastable AFE states are not energy accessibleby experimental E, the CCPS lattice spontaneously responds to theexternal electric field by readjusting the positions of Cu+ ions withinthe vdW cages of CrS6 and P2S6 octahedrons. Our DFT calculationsreveal that, the polarizations are predominated by layer-dependentrearrangements of the anti-parallel Cu+ ions, leading the formationof a heterogenous FiE state. As schemed in the second column ofFig. 4a, due to the unique vdW depolarization effect, anti-parallelCu+ ions in different MLs will move toward the Cr atomic planeunder E, with Δd proportional to the energy of the ML-indexedmetastable AFE states. The formation of an E-dependent FiE state isconfirmed by SHG experiments, with tunable Vb applied to theheterostructures of FLG/BN/CCPS/BN/FLG. As shown in Fig. 4c, bycooling down a 6L device, the T-dependent SHG intensity decreasesmonotonically and eventually stabilizes at about 140 K, whichmanifests the AFE transition of thin-layer CCPS when Cu+ ions lockinto the stripe-AFE positions. We also conduct polarization-dependent SHG measurements at 300 K and 77 K, respectively,with different Vb setpoints. Strikingly, polarization-SHG at 77 Kunveils a fingerprinting six-fold to two-fold symmetry transitionwhen thin-layer CCPS under E is driven from the AFE state into theFiE phase (Fig. 4d). In stark contrast, the same CCPS device at 300 Kconsistently shows a nearly six-fold symmetry for different Vb(Supplementary Fig. 14c). Similar results are well reproduced onmultiple devices, and also repeated on an FLG/CCPS/SiO2/Si struc-ture, in which the heavily doped silicon substrate was used as thebias electrode. The Vb-dependent SHG intensity reflects continuousinversion symmetry breaking induced by E. As illustrated in theinsets of Fig. 4d, for the stripe-AFE ground state, there exists aninversion center between two neighboring monolayers, while theVb-driven heterogeneous FiE state has reduced inversion symmetrydue to layer-dependent rearrangements of anti-parallel Cu+ ionswithin the FM vdW cages.Fig. 4 | Ferrielectricity controlledmagnetoelectric coupling in thin-layer CCPS.a Latticemodel representationof the stripe-AFE state (first column), the E-enforcedheterogenous FiE phase (second column), and energy unfavorable metastable AFEstates (third column; grouped in orange dashed lines) for BL and trilayer CCPS,respectively. Here, the vdW monolayer cages are represented by blue/orange rec-tangles, standing for energy accessible and inaccessible respectively. Within thecages, the ground-state (polarized) Cu+ ions are depicted byblue (red) circles. Here,2L-mAFE-I and 2L-mAFE-II correspond to the flipping of an anti-parallel Cu-ionstripe of within the bottom and top ML, respectively. b Black curve: Energypotential diagram vs Δd (anti-parallel Cu+ displacement) of BL CCPS from the AFEground state to 2L-mAFE-I. Red curve:ΔEMEConΔd, whichmonotonically decreaseswhenCu+ ionsmove towards the Cr atomic plane. c SHG intensity vs T, showing theAFE transition and the frozen ofCu+ ions at around 140 K.d Polarization-dependentSHG at 77K under different Vb. The distinctive six-fold to two-fold transition in thepolarization-SHG patterns proves that the non-linear optical anisotropy of thin-layer CCPS associated with inversion symmetry breaking can be continuouslytuned by Vb. e Illustration of the FiE-interlockedMECmechanism in CCPS, in whichlattice distortions introduced by Cu+ ion displacements efficiently lower theinterlayer superexchange interaction.Article https://doi.org/10.1038/s41467-024-47373-7Nature Communications |         (2024) 15:3029 5Most unusually, the formation of E-dependent heterogenous FiEstate is accompanied by a significant reduction in themagnetic energydifference between the interlayer FM0 and AFM configurations, whichmanifests the change of Vb-tunable MEC energy ΔEMEC. Using theillustrative example of BL CCPS, the DFT calculations deduce that theFiE state has monotonically decreased ΔEMEC(E), defined byΔEMECðEÞ= ðEFM0 ðEÞ � EAFMðEÞÞ � ðEFM0 ð0Þ � EAFMð0ÞÞ, when E drivesanti-parallel positioned Cu+ ions toward the energy saddle point of theCr atomic plane (see Supplementary Note 7 for the detailed formula-tion of ΔEMEC(E)). As illustrated in Fig. 4e, for thin-layer CCPS, theinterlayer superexchange that determines J⊥ is intermediated by theinterfacial majority spin p-orbitals of S atoms, which are antiparallel tothe Cr magnetic moments due to the formation of Cr-S valencebonds57. The stripe-AFE ground state corresponds to a relative sym-metric CrS6 octahedron with minimal lattice distortions (Fig. 4e andDFT calculations in Supplementary Table 1). In response to externalelectric field, anti-parallel Cu+ ions move within the FM vdW cages,enforcing substantially lattice distortions for the CrS6 octahedrons bystretching the average Cr-S bond lengths. The resulting attenuation ofthe Cr-S bonding energy effectively decreases the netmajority spins ofS atoms, which changes from 0.135 electrons per S for d =0Å (theequilibrium stripe-AFE position) to 0.092 electrons per S forΔd =0.44Å, and thus, efficiently lowers the interlayer AFM coupling.As summarized in Fig. 4b, the collective movements of anti-parallel Cu+ ions within the FM vdW cages reach the maximum ΔEMECof −0.37meV for the saddle point position Δd = 1.69Å. It should benoted that ΔEMEC(E) has a maximized change rate within the E-polar-izable regime, which is highlighted by the blue color (E < 1.5 V nm−1; seeSupplementary Note 6 for electric breakdown test). In this polarizationregime, ΔEMEC rapidly reaches −0.15meV for a minute Δd of 0.29Å,which is qualatatively in agreement with the experimental ΔEMEC esti-mated by H ±sat. As a direct consequence of this ultra Vb-tunable ΔEMEC,the Vb-dependent heterogenous FiE state has drastically reduced J⊥,which changes by nearly 50%with amoderate E of 0.3V nm−1 as shownin Fig. 3c. After elucidating the Vb-tunable MEC mechanism, the finalkey element is to understand the Vb-polarity asymmetry in the MECcoupling, which is attributed to a build-in potential introduced by theasymmetric work functions of top and bottom FLG electrodes. UsingFowler-Nordheim (FN) tunneling in complementary with scanningKelvin probe microscope (SKPM), we determined that the asymmetricwork functions can readily reach 0.5 eV due to the Dirac electronicstructure of graphene and the interfacial doping byAFECCPS flakes55,56(see Supplementary Note 4 and Note 8 for detailed analyses using FNtunneling and SKPM). As shown in Supplementary Fig. 7a, the exis-tence of a build-in potential makes the two otherwise Vb-polarityequivalent FiE states asymmetric in energy, in analogy to GFeFETsunder a constant electrostatic doping55,56.DiscussionThe striking energy difference between the 2L-mAFE-I and 2L-mAFE-IIstates of BL CCPS provide an amazing example on how the well-established depolarization mechanism may become peculiarly differ-ent at the 2D limit in a vdW multiferroic. The discovered FiE inter-locked widely-tunable MEC mechanism also suggests the greatpotential of realizing highly-efficient electric control of 2D magnetismin emergent vdW multiferroics. Equally important, the multiferroicityin thin-layer CCPS opens the possibility for the investigation of 2Dquantum phase transitions, and the fabrication of complex vdW het-erostructure for the realization of two-phase multiferroic for novelelectronically controlled spintronics and Josephson junctions.MethodsCrystal growthBulk CCPS single crystals were synthesized by the standard self-fluxmethod. High purity Cu (99.999%), Cr (99.996%), P (99.999%) and S(99.5%) powders were mixed stoichiometrically and loaded into aquartz tube in an argon-filled glove box. The quartz tube was evac-uated to a pressure of 10−2 Pa before flame-sealing and placing into atube furnace. For crystal growth, the temperature was ramped fromroom-temperature to 973 K, which was held constant for 14 days byPID control. After that, the furnace was naturally cooled to roomtemperature.Tunneling junction fabricationsCCPS MTJs were fabricated in an argon-filled glove box. Differentthicknesses of CCPS flakes were exfoliated from single-crystal seeds,identified by optical contrast and further cross checked by atomicforce microscopy (Park system NX10) and Raman spectroscopy. Toobtain aMTJ shown in Fig. 1e, thin-layer CCPS flakes were exfoliated onSiO2/Si substrates, followed by two PC/PDMS based dry-transfersteps58 to prepare the top and bottom few-layer graphene (FLG) elec-trodes. In the final step, the MTJ is covered by a thin-film flake ofhexagonal boron nitride to protect the sample from ambient expo-sure. Top and bottom FLG electrodes were contacted by 5/50nmCr/Au electrode using standard electron beam lithography (EBL)technique followed by thermal evaporation.Electrical measurementsElectrical measurements were performed in an Oxford-14 T cryostat.The system is equipped with a sample rotator for applying magneticfield either in-plane or out-of-plane. The I −V curve were measuredusing a Keithley 2400 and zero-bias-resistance were measure by usingstandard lock-in method (with excitation frequency < 20Hz).SHG measurementsThe SHG measurements are conducted in a liquid nitrogen cooledoptical cryostat with a femtosecond laser (Rainbow OEM, NPI Lasers).The wavelength is centered at either 1064 nm or 1560nm, bothyielding consistent SHG results. The laser beam is focused on thesamples with a near-infrared radiation (NIR) objective (MplanApo NIR50x, OptoSigma) and collected in a reflective geometry, wherein theSHG signal is separated with a dichroic mirror. The signal is detectedwith a fiber coupled spectrometer (UHTS 600 VIS, Witec). Forpolarization-dependent SHG measurements, the polarization anglebetween the incident laser and the SHG signal is set by two polarizers,and the in-plane polarization direction is rotated by a dual band half-wave plate before the optical objective. The optic setup is schemati-cally illustrated in Supplementary Fig. 14d. The signals are confirmedto be originating in the CCPS samples by comparingwith the referencespectra from BN, FLG electrodes and SiO2 areas on the substrates.DFT calculationsThe first-principles DFT calculations were performed using the Viennaab initio simulation package (VASP) with the choice of projector aug-mented waves (PAW) basis set59–61. The Perdew-Burke-Ernzerh (PBE)functional has been employed to treat the exchange and correlationfunctional62. The energy cut off for the plane-wave basis was set to600 eV. All the structures were fully relaxed until the force on everyatom and energy were converged to 0.01 eÅ−1 and 1 × 10−7 eV, respec-tively. A Γ-centred Monkhorst-Pack 12 × 8 × 1 k-mesh was used fork-point sampling63. To model the 2D films, a vacuum layer of at least20Å was included in the c-axis direction. The Hubbard correction ofU = 3 eV is applied to the 3d orbitals of Cr64. The vdW corrections wereincluded by the DFT-D3 method65. The kinetic pathways were calcu-lated by the climbing image nudged elastic band (CINEB) method66.Layer-dependent CAFM saturation field HsatFor in-plane H∥, the magnetization of each ML (Ms) lies within the vdWplane, which can be described by a 1D spin-chain model24. By labelingthemagnetization of the ith layer asMi =Msðcosθi, sinθi,0Þ, whereθi isArticle https://doi.org/10.1038/s41467-024-47373-7Nature Communications |         (2024) 15:3029 6the in-plane CAFM canting angle, the magnetic coupling energy of anN-layer CCPS is:UNðθ1, . . . ,θN;HÞ= JXN�1i = 1cosðθi + 1 � θiÞ � μ0MsHXNi = 1cosðθiÞ: ð1ÞThe energy favorable CAFM state for a given H can be obtained bysolving the lowest-energy solutions of UN,∂UN∂θi=0 8i= 1, . . . , N : ð2ÞFor a solution set of θsi = 1, . . . , N� �satisfying the energy minimumrequirement of UN, the matrix of the second derivatives of UN musthave all positive eigenvalues:det ðAsÞ>0 with Asij =∂2UN∂θi∂θj�����θi =θsi: ð3ÞFor the saturation FM0 state with θ1 = θ2 =… = θN = 0, the correspond-ing matrix AFM0has a much simplified tridiagonal form of,AFM 0=H � HJ2HJ2 0HJ2 H � HJHJ2HJ2 . . . . . .. . . . . .HJ2HJ2 H � HJHJ20 HJ2 H � HJ20BBBBBBBBBBBB@1CCCCCCCCCCCCA, ð4Þwhere HJ = 2J/(μ0Ms) represents the energy induced by the interlayerexchange coupling J (J⊥). For AFM0, the (N − 1) eigenvalues can beobtained by the standard approach for tridiagonal matrices, whichreads:λFM0k =H � HJ +HJ coskπN� �with k = 1, . . . , N � 1: ð5ÞTo get a positive matrix determinant, Hsat must make the smallesteigenvalue for k =N − 1 vanish:Hsat � HJ +HJ cosN� 1Nπ� �=0, ð6Þwhich yield the analytical result of Hsat byHsat =HJ +HJ cosπN� �=2HJcos2 π2N� �: ð7ÞEq. (7) explains the layer-dependent CAFM saturation field Hsat, whichincreases substantially from 3.2 T for BL CCPS to 5.7 T for quintuple-layer CCPS. Using experimental data, the interlayer exchange couplingJ can also be extracted by fitting the layer-dependent Hsat to Eq. (7).Data availabilityThe data that support the findings of this study are provided in theSource Data file. Source data are provided with this paper.References1. Eerenstein, W., Mathur, N. D. & Scott, J. F. Multiferroic and mag-netoelectric materials. Nature 442, 759–765 (2006).2. Wang, J. et al. Epitaxial BiFeO3 multiferroic thin film hetero-structures. Science 299, 1719–1722 (2003).3. Kimura, T. et al. Magnetic control of ferroelectric polarization. Nat-ure 426, 55–58 (2003).4. Khomskii, D. Classifying multiferroics: mechanisms and effects.Physics 2, 20 (2009).5. Cheong, S.-W. & Mostovoy, M. Multiferroics: a magnetic twist forferroelectricity. Nat. Mater. 6, 13–20 (2007).6. Huang, B. et al. Layer-dependent ferromagnetism in a van derWaals crystal down to the monolayer limit. Nature 546, 270–273(2017).7. Gong, C. et al. Discovery of intrinsic ferromagnetism in two-dimensional van der Waals crystals. Nature 546, 265 (2017).8. Chang, K. et al. Discovery of robust in-plane ferroelectricity inatomic-thick SnTe. Science 353, 274–278 (2016).9. Song, Q. et al. Evidence for a single-layer van der Waals multi-ferroic. Nature 602, 601–605 (2022).10. Gou, J. et al. Two-dimensional ferroelectricity in a single-elementbismuth monolayer. Nature 617, 67–72 (2023).11. Wang, Z. et al. Electric-field control of magnetism in a few-layeredvan derWaals ferromagnetic semiconductor.Nat. Nanotechnol. 13,554–559 (2018).12. Wang, Y. et al. Strain-sensitive magnetization reversal of a van derWaals magnet. Adv. Mater. 32, 2004533 (2020).13. Deng, Y. et al. Gate-tunable room-temperature ferromagnetism intwo-dimensional Fe3GeTe2. Nature 563, 94–99 (2018).14. Zhang,M. J. et al. Spin-lattice coupledmetamagnetism in frustratedvan der Waals magnet CrOCl. Small 19, 2300964 (2023).15. Wang, J. et al. Giant magnetoelectric effect at the graphone/fer-roelectric interface. Sci. Rep. 8, 12448 (2018).16. Shang, J. et al. Robust magnetoelectric effect in the decoratedgraphene/In2Se3 heterostructure. Adv. Mater. Int. 13,3033–3039 (2021).17. Wang, Z. et al. Very large tunneling magnetoresistance in layeredmagnetic semiconductor CrI3. Nat. Commun. 9, 2516 (2018).18. Kim, H. H. et al. One million percent tunnel magnetoresistance in amagnetic van der Waals heterostructure. Nano Lett. 18,4885–4890 (2018).19. Ghazaryan, D. et al. Magnon-assisted tunnelling in van der Waalsheterostructures based on CrBr3. Nat. Electron. 1, 344–349 (2018).20. Fei, Z. et al. Two-dimensional itinerant Ising ferromagnetism inatomically thin Fe3GeTe2. Nat. Mater. 17, 778–782 (2018).21. Kim, K. et al. Suppression of magnetic ordering in XXZ-type anti-ferromagnetic monolayer NiPS3. Nat. Commun. 10, 345 (2019).22. Cai, X. et al. Atomically thin CrCl3: an in-plane layered anti-ferromagnetic insulator. Nano Lett. 19, 3993–3998 (2019).23. Klein, D. R. et al. Enhancement of interlayer exchange in an ultrathintwo-dimensional magnet. Nat. Phys. 15, 1255–1260 (2019).24. Wang, Z. et al. Determining the phase diagram of atomically thinlayered antiferromagnet CrCl3. Nat. Nanotechnol. 14,1116–1122 (2019).25. Long, G. et al. Persistence of magnetism in atomically thin MnPS3crystals. Nano Lett. 20, 2452–2459 (2020).26. Telford, E. J. et al. Layered antiferromagnetism induces largenegative magnetoresistance in the van der Waals semiconductorCrSBr. Adv. Mater. 32, 2003240 (2020).27. Peng, Y. et al. Magnetic structure and metamagnetic transitions inthe van der Waals antiferromagnet CrPS4. Adv. Mater. 32,202001200 (2020).28. Kim, H. H. et al. Magneto-memristive switching in a 2D layer anti-ferromagnet. Adv. Mater. 32, 1905433 (2020).29. Liu, F. et al. Room-temperature ferroelectricity in CuInP2S6 ultrathinflakes. Nat. Commun. 7, 12357 (2016).30. Xiao, J. et al. Intrinsic two-dimensional ferroelectricity with dipolelocking. Phys. Rev. Lett. 120, 227601 (2018).31. Wan, S. et al. Non-volatile ferroelectricmemory effect in ultrathin α-In2Se3. Adv. Funct. Mater. 29, 1808606 (2019).Article https://doi.org/10.1038/s41467-024-47373-7Nature Communications |         (2024) 15:3029 732. Fei, Z. et al. Ferroelectric switching of a two-dimensional metal.Nature 560, 336–339 (2018).33. Yasuda, K., Wang, X., Watanabe, K., Taniguchi, T. & Jarillo-Herrero,P. Stacking-engineered ferroelectricity in bilayer boron nitride.Science 372, 1458–1462 (2021).34. Stern,M. V. et al. Interfacial ferroelectricity byvanderWaals sliding.Science 372, 1462–1466 (2021).35. Fei, R., Kang, W. & Yang, L. Ferroelectricity and phase transitions inmonolayer group-IV monochalcogenides. Phys. Rev. Lett. 117,097601 (2016).36. Wang, H. & Qian, X. Two-dimensional multiferroics in monolayergroup IV monochalcogenides. 2D Mater. 4, 015042 (2017).37. Seixas, L., Rodin, A. S., Carvalho, A. &Castro Neto, A. H.Multiferroictwo-dimensional materials. Phys. Rev. Lett. 116, 206803 (2016).38. Susner, M. A., Chyasnavichyus, M., McGuire, M. A., Ganesh, P. &Maksymovych, P. Metal thio- and selenophosphates as multi-functional van der Waals layered materials. Adv. Mater. 29,1602852 (2017).39. Lee, J.-U. et al. Ising-type magnetic ordering in atomically thinFePS3. Nano Lett. 16, 7433–7438 (2016).40. Colombet, P., Leblanc, A., Danot, M. & Rouxel, J. Structural aspectsand magnetic properties of the lamellar compoundCu0.50Cr0.50PS3. J. Solid State Chem. 41, 174–184 (1982).41. Maisonneuve, V., Payen, C. & Cajipe, V. B. On CuCrP2S6 : copperdisorder, stacking distortions, andmagnetic ordering. J. Solid StateChem. 116, 208–210 (1995).42. Kleemann, W., Shvartsman, V. V., Borisov, P., Banys, J. & Vyso-chanskii, Y.M.Magnetic and polar phases anddynamical clusteringin multiferroic layered solid solutions CuCr1−xInxP2S6. Phys. Rev. B84, 094411 (2011).43. Dziaugys, A. et al. Phase diagramofmixedCu(InxCr1−x)P2S6 crystals.Phys. Rev. B 85, 134105 (2012).44. Lai, Y. et al. Two-dimensional ferromagnetism and driven ferroe-lectricity in van der Waals CuCrP2S6. Nanoscale 11, 5163–5170(2019).45. Selter, S. et al. Crystal growth, exfoliation, andmagnetic propertiesof quaternary quasi-two-dimensional CuCrP2S6. Phys. Rev. Mater. 7,033402 (2023).46. Wang, X. et al. Electrical andmagnetic anisotropies in vanderWaalsmultiferroic CuCrP2S6. Nat. Commun. 14, 840 (2023).47. Ma, R.-R. et al. Nanoscale mapping of Cu-ion transport in van derWaals layered CuCrP2S6. Adv. Mater. Int. 9, 2101769 (2022).48. Park, C. B. et al. Observation of spin-induced ferroelectricity in alayered van der Waals antiferromagnet CuCrP2S6. Adv. Electron.Mater. 8, 2101072 (2022).49. Cho, K. et al. Tunable ferroelectricity in van der Waals layeredantiferroelectric CuCrP2S6. Adv. Funct. Mater. 32, 2204214 (2022).50. Hua, C. et al. Strong coupled magnetic and electric ordering inmonolayer of metal thio(seleno)phosphates. Chin. Phys. Lett. 38,077501 (2021).51. Io, W. F. et al. Direct observation of intrinsic room-temperatureferroelectricity in 2D layered CuCrP2S6. Nat. Commun. 14,7304 (2023).52. Worledge, D. C. & Geballe, T. H. Magnetoresistive double spin filtertunnel junction. J. Appl. Phys. 88, 5277–5279 (2000).53. Klein, D. R. et al. Probingmagnetism in 2D van derWaals crystallineinsulators via electron tunneling. Science 360, 1218–1222 (2018).54. Song, T. et al. Giant tunneling magnetoresistance in spin-filter vander Waals heterostructures. Science 360, 1214–1218 (2018).55. Zheng, Y. et al. Gate-controlled nonvolatile graphene-ferroelectricmemory. Appl. Phys. Lett. 94, 163505 (2009).56. Zheng, Y. et al. Graphene field-effect transistors with ferroelectricgating. Phys. Rev. Lett. 105, 166602 (2010).57. Wang, C. et al. Bethe-slater-curve-like behavior and interlayer spin-exchange coupling mechanisms in two-dimensional magneticbilayers. Phys. Rev. B 102, 020402 (2020).58. Zomer, P. J., Guimar aes,M., Brant, J. C., Tombros, N. & VanWees, B.J. Fast pick up technique for high quality heterostructures of bilayergraphene and hexagonal boron nitride. Appl. Phys. Lett. 105,013101 (2014).59. Kresse, G. & Furthmüller, J. Efficient iterative schemes for ab initiototal-energy calculations using a plane-wave basis set. Phys. Rev. B54, 11169 (1996).60. Kresse, G. & Furthmüller, J. Efficiency of ab-initio total energy cal-culations for metals and semiconductors using a plane-wave basisset. Comput. Mater. Sci. 6, 15–50 (1996).61. Blöchl, P. E. Projector augmented-wave method. Phys. Rev. B 50,17953 (1994).62. Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized gradientapproximation made simple. Phys. Rev. Lett. 77, 3865 (1996).63. Monkhorst, H. J. & Pack, J. D. Special points for Brillouin-zoneintegrations. Phys. Rev. B 13, 5188 (1976).64. Dudarev, S. L., Botton, G. A., Savrasov, S. Y., Humphreys, C. J. &Sutton, A. P. Electron-energy-loss spectra and the structural stabi-lity of nickel oxide: an LSDA+U study. Phys. Rev. B 57, 1505 (1998).65. Grimme, S., Ehrlich, S. & Goerigk, L. Effect of the damping functionin dispersion corrected density functional theory. J. Comput. Chem.32, 1456 (2011).66. Henkelman, G., Uberuaga, B. P. & Jónsson, H. A climbing imagenudged elastic band method for finding saddle points and mini-mum energy paths. J. Chem. Phys. 113, 9901–9904 (2000).AcknowledgementsThis work is supported by the National Key R&D Program of theMOST ofChina (Grant Nos. 2023YFA1406302, and 2019YFA0308602), theNational Science Foundation of China (Grant Nos. 12374194, 11574264,12241401 and 12174336), and the Zhejiang Provincial Natural ScienceFoundation (D19A040001). Y.Z. acknowledges support from the Userswith Excellence Project of Hefei Science Center CAS, 2021HSC-UE007.This work made use of the resources of the Center of Electron Micro-scopy of Zhejiang University.Author contributionsY.Z. initiated and supervised the project. Q.F.H. and M.J.Z. synthesizedand characterized CCPS crystals. Q.F.H., M.J.Z. and Y.W. fabricatedCCPS devices and carried out all the measurements. S.J.D. and C.H.J.did the ADF-STEM experiments. Y.Q.H., C.Q.H. and Y.H.L. did the DFTcalculations. T.T. and K.W. synthesized the hBN crystals. X.F.X., J.B.Y.,S.J.Y. and D.W.W. discussed the results and commented on the manu-script. Q.F.H., Y.Q.H., Y.W., L.J.L. and Y.Z. analysed the data and wrotethe paper with inputs from all authors.Competing interestsThe authors declare no competing interests.Additional informationSupplementary information The online version containssupplementary material available athttps://doi.org/10.1038/s41467-024-47373-7.Correspondence and requests for materials should be addressed toLinjun Li, Yunhao Lu, Chuanhong Jin or Yi Zheng.Peer review information Nature Communications thanks BingchengLuo, Soonyong Park and Hailin Peng for their contribution to the peerreview of this work. A peer review file is available.Article https://doi.org/10.1038/s41467-024-47373-7Nature Communications |         (2024) 15:3029 8https://doi.org/10.1038/s41467-024-47373-7Reprints and permissions information is available athttp://www.nature.com/reprintsPublisher’s note Springer Nature remains neutral with regard to jur-isdictional claims in published maps and institutional affiliations.Open Access This article is licensed under a Creative CommonsAttribution 4.0 International License, which permits use, sharing,adaptation, distribution and reproduction in any medium or format, aslong as you give appropriate credit to the original author(s) and thesource, provide a link to the Creative Commons licence, and indicate ifchanges were made. The images or other third party material in thisarticle are included in the article’s Creative Commons licence, unlessindicated otherwise in a credit line to the material. If material is notincluded in the article’s Creative Commons licence and your intendeduse is not permitted by statutory regulation or exceeds the permitteduse, you will need to obtain permission directly from the copyrightholder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.© The Author(s) 20241School of Physics, and State Key Laboratory of Silicon and Advanced Semiconductor Materials, Zhejiang University, Hangzhou 310027, China. 2ZhejiangProvince Key Laboratory of Quantum Technology and Device, School of Physics, Zhejiang University, Hangzhou 310027, China. 3State Key Laboratory ofSilicon and Advanced Semiconductor Materials, Zhejiang University, Hangzhou 310027, China. 4State Key Laboratory of Extreme Photonics and Instru-mentation, College of Optical Science and Engineering, Zhejiang University, Hangzhou 310027, China. 5School of Physics Science and Engineering, TongjiUniversity, Shanghai 200092, China. 6State Key Laboratory for Mesoscopic Physics, School of Physics, Peking University, Beijing 100871, China. 7School ofPhysics and Technology, Wuhan University, Wuhan 430072, China. 8National Institute for Materials Science, 1-1 Namiki, Tsukuba 305-0044, Japan. 9Theseauthors contributed equally: Qifeng Hu, Yuqiang Huang, Yang Wang, Sujuan Ding, Minjie Zhang. e-mail: lilinjun@zju.edu.cn; luyh@zju.edu.cn;chhjin@zju.edu.cn; phyzhengyi@zju.edu.cnArticle https://doi.org/10.1038/s41467-024-47373-7Nature Communications |         (2024) 15:3029 9http://www.nature.com/reprintshttp://creativecommons.org/licenses/by/4.0/http://creativecommons.org/licenses/by/4.0/mailto:lilinjun@zju.edu.cnmailto:luyh@zju.edu.cnmailto:chhjin@zju.edu.cnmailto:phyzhengyi@zju.edu.cn Ferrielectricity controlled widely-tunable magnetoelectric coupling in van der Waals multiferroics Results Discussion Methods Crystal�growth Tunneling junction fabrications Electrical measurements SHG measurements DFT calculations Layer-dependent CAFM saturation field Hsat Data availability References Acknowledgements Author contributions Competing interests Additional information