\section*{Abstract}
{A promising approach to the next generation of low-power, functional, and energy-efficient electronics relies on novel materials with coupled magnetic and electric degrees of freedom. In particular, stripy antiferromagnets often exhibit broken crystal and magnetic symmetries, which may bring about the magnetoelectric (ME) effect and enable the manipulation of intriguing properties and functionalities by electrical means. The demand for expanding the boundaries of data storage and processing technologies has led to the development of spintronics toward two-dimensional (2D) platforms. This work reports the ME effect in the 2D stripy antiferromagnetic insulator CrOCl down to a single layer. By measuring the tunneling resistance of CrOCl on the parameter space of temperature, magnetic field, and applied voltage, we verified the ME coupling down to the 2D limit and unraveled its mechanism. Utilizing the multi-stable states and ME coupling at magnetic phase transitions, we realize multi-state data storage in the tunneling devices. Our work not only advances the fundamental understanding of spin-charge coupling, but also demonstrates the great potential of 2D antiferromagnetic materials to deliver devices and circuits beyond the traditional binary operations.}
\end{quotation}
}]


\newpage 
\clearpage



\noindent 
\section*{Introduction}
\vspace{-2mm}

The field of spintronics concerns the processing of digital information, where an external stimulus, preferably an electrical stimulus, is applied to control the spin order in a magnetic system that acts as a "0" or "1" digital bit. As a more common magnetic ground state than ferromagnetism, antiferromagnetism has received increasing attention in recent years due to its promising prospect for spintronic devices, such as their robustness against external perturbations, no stray fields, and ultrafast dynamics \cite{jungwirth2016antiferromagnetic}. More intriguingly, antiferromagnets often manifest complicated spin configurations and phase transitions, opening up the possibility to implement new data storage logic superior to conventional binary algorithms. However, the negligible net magnetization also makes it difficult to read out or electrically manipulate the antiferromagnetic order, which is desired for information technology. As a result, formidable efforts have been devoted to switching antiferromagnets through spin torques from exchange bias\cite{fukami2016magnetization}, spin-orbit torque\cite{tsai2020electrical} and spin-galvanic effect\cite{wadley2016electrical}, etc.\par

Apart from the above-mentioned "external" approaches, the "internal" approach to antiferromagnetic spintronics rests upon materials with coupled degrees of freedom \cite{matsukura2015control,eerenstein2006multiferroic,spaldin2019advances} such as magnetic moment, electric polarization and strain, which are often referred to as the magnetoelastic or magnetoelectric materials\cite{yang2015bifeo3}. However, as the transition metal $d$ electrons are supposed to repel the tendency for off-center ferroelectric distortion\cite{hill2000there}, the coexistence of magnetic moment and electric polarization is hard to achieve. The coupled energy terms, as well as the electric polarization, require materials with low crystal and magnetic symmetry. Improper magnetic ferroelectrics\cite{cheong2007multiferroics,tokura2014multiferroics} where the ferroelectricity originated from spin-order-driven inversion symmetry breaking are considered to be an ideal platform to realize the mutual clamping of the order parameters of ferroelectricity and antiferromagnetism\cite{fiebig2002observation,lottermoser2004magnetic,lorenz2004field}. Based on the spin frustration and spin-orbit coupling theories, three spin-structure-induced electric polarization mechanisms have been established, namely the exchange striction model\cite{sergienko2006ferroelectricity}, the inverse Dzyaloshinskii-Moriya interaction model\cite{katsura2005spin} and the \textit{p-d} hybridization model\cite{jia2006bond}.

To further explore ME-based fundamental physics and develop practical applications for information devices, it is inevitable to extend the research to the two-dimensional (2D) limit\cite{lin2019two,sierra2021van,huang2018electrical,jiang2018electric,sharpe2019emergent,he2020giant}. Recently, multiferroicity resulting from inverse Dzyaloshinskii-Moriya interaction\cite{ju2021possible,song2022evidence} and \textit{p-d} hybridization\cite{park2022observation} have been evidenced in van der Waals materials. Nevertheless, the direct coupling between magnetic moment and electric polarization in the 2D limit has not been reported. The air-stable van der Waals insulator CrOCl, with a stripy antiferromagnetic ground state and thus both broken rotational symmetry and translation symmetry perpendicular to the stripes, is a potential candidate to exhibit the intrinsic 2D ME effect. We have verified the magnetic phase transitions and magnetoelastic coupling effect of CrOCl in the 2D limit\cite{gu2022magnetic}. Additionally, exotic quantum Hall effects\cite{wang2021flavoured} and insulator phases\cite{yang2021realization} were observed in the graphene/CrOCl heterostructures. In this work, we report the intrinsic ME effect in CrOCl down to the 2D limit. By means of dielectric measurements combined with first-principles calculations, we find that the wave symmetry of stripy antiferromagnetism induces antiferroelectric dipoles below the Néel temperature. The antiferroelectric dipoles persist to the 2D limit, as evidenced by the $I-V$ hysteresis loops and magnetoresistance of CrOCl tunneling devices. We also realize direct manipulation of resistance states through electric and magnetic excitations, demonstrating a data storage device with continuously adjustable outputs. The metastable states in the first-order phase transition resulting from the ME coupling break through the barriers of write-in and read-out information in antiferromagnets and hopefully can open an effective route for multi-state memories with high storage density, nonvolatility, low energy consumption, and small memory cell.

\vspace{-2mm}
\section*{Results and Discussion}
\vspace{-1mm}

\begin{figure}[ht!]
\includegraphics[width=0.48\textwidth]{Figure1.jpg}
\caption{\textbf{Magnetic and dielectric properties of bulk CrOCl.} 
\textbf{(a),} Crystal structure of CrOCl above the Néel temperature with the $Pmmn$ space group in a unit cell. \textbf{(b),} Crystal structure of CrOCl below the Néel temperature with the $P2_1/m$ space group in a unit cell. The structural parameter $\alpha$ that defines the crystal symmetry is labeled in the unit cell. The vectors on each atom represent the atomic displacements (magnified 100 times), while the dashed boxes and open arrows represent the calculated antiferroelectric polarizations of different Cl-Cr-O chains. \textbf{(c),} Temperature dependence of relative permittivity $\epsilon_r$ and magnetic moment of bulk CrOCl crystal. \textbf{(d),} $\Delta\epsilon_r/\epsilon_r$(0 T) \textit{versus} out-of-plane magnetic field and the corresponding $M-H$ curve at 2 K.
}
\label{F1}
\end{figure}

\begin{figure*}[ht!]
\centering 
\includegraphics[width=1\textwidth]{Figure2.jpg}
\caption{\textbf{$I-V$ hysteresis of CrOCl tunneling devices.} 
\textbf{(a),} Illustration of the tunneling device. \textbf{(b),} $I-V$ curves of a few-layer CrOCl ($\sim$10 nm) tunneling device at 20 K (above the Néel temperature) and 2 K (below the Néel temperature). \textbf{(c),} 2D map of current polarization $\rho$ as a function of temperature and bias voltage. The green dashed line marks the Néel temperature of exfoliated CrOCl. \textbf{(d),} $\rho$ \textit{versus} temperature at different voltages extracted from (b). \textbf{(e),} Resistance \textit{versus} temperature of CrOCl tunneling device at 5 V in different sweeping processes. \textbf{(f),} Tunneling resistance \textit{versus} out-of-plane magnetic field at $-$5 V in different sweeping processes. \textbf{(g),} Reproducible manipulation of the resistance states at $-$0.8 V by alternately changing the sweeping direction. Data in (a-f) were obtained in the $\sim$10 nm CrOCl tunneling device, while (g) was obtained in a single-layer CrOCl tunneling device in parallel with a 10 M$\Omega$ protection resistor.
}
\label{F2}
\end{figure*}

\begin{figure*}[ht!]
\centering 
\includegraphics[width=0.85\textwidth]{Figure3.jpg}
\caption{\textbf{Operation principles of multi-state data storage.} 
\textbf{(a),}$I-B$ curve of the CrOCl tunneling device at 2 K. The pink curve shows the differential current $\delta I=I(B)|_{\rm{up}}-I(B)|_{\rm{down}}$. The critical field where $\delta I$ reaches a maximum ($B_0$) and the hysteresis closes ($B_s$) are annotated by black arrows. \textbf{(b),} Tunneling current over time after several electric excitations. The dashed lines represent the current values at 3.33 T in $B_{\rm{down}}$ and $B_{\rm{up}}$ sweeping cycles extracted from (a). \textbf{(c),} $I-B$ curves with different magnetic sweep ranges. \textbf{(d),} 2D map of the tunneling current after designed electric and magnetic excitations. \textbf{(e,f),} Tunneling current (e) after alternating electric and magnetic excitations. The corresponding magnetic and electric fields \textit{versus} time are plotted in (f).
}
\label{F3}
\end{figure*}

To understand the intrinsic ME effect, we start with the energy coupling terms for the ground state of CrOCl. As reported by several neutron scattering results \cite{norlund1975preparation,reuvekamp2014investigation,zhang2014temperature} and confirmed by DFT calculations\cite{gu2022magnetic}, CrOCl transforms into the $\uparrow \uparrow \downarrow \downarrow$ stripy AFM ground state at $\sim$14 K. The magnetic frustration drives the crystal from the orthorhombic space group $Pmmn$ (Fig.\ref{F1}a) to the monoclinic space group $P2_1/m$ (Fig. \ref{F1}b) \cite{angelkort2009observation,reuvekamp2014investigation,zhang2014temperature}. Despite the monoclinic distortion, the space inversion operation $P$ preserves in the $\uparrow \uparrow \downarrow \downarrow$ ground state due to the special zigzag atom-chain structure of CrOCl. The two sublayers of a single van der Waals layer are spatially inverse to each other, thus preventing net electric polarization from appearing. However, based on the symmetry operations of a wave order, we can obtain the free energy in the form\cite{balatsky2010induction,zachar1998landau}:

\vspace{-2mm}
\begin{equation}
\begin{aligned}
 \mathcal{F}=& \frac{1}{2} r_{a}|A_\mathbf{k}|^{2}+\frac{1}{2} r_{s}|\mathbf{S_q}|^{2}\\
 &+\lambda_{1}\left[(\mathbf{S_q} \cdot \mathbf{S_q}) A_\mathbf{k}^{*}+\text { c.c. }\right]+2 \lambda_{2}|\mathbf{S_q}|^{2}|A_\mathbf{k}|^{2}
\end{aligned}
\label{EQ1}
\end{equation}

\noindent where $\mathbf{S_q}$ is the amplitude of spin density modulation with the wave vector of $\mathbf{q}$ and $A_\mathbf{k}$ is the corresponding lattice modulation with the wave vector of $\mathbf{k}$. $r_a$, $r_s$, $\lambda_1$ and $\lambda_2$ are constant parameters of the material. Minimizing the free energy, we obtain the lattice modulation in a collinear spin configuration:

\vspace{-2mm}
\begin{equation}
A_\mathbf{k}=\frac{\lambda_1}{r_a} \mathbf{S_q} \cdot \mathbf{S_q}+\text {c.c.}
\label{EQ2}
\end{equation}

\noindent The lattice modulation, consequently, can only be present with $\mathbf{k}=\pm2\mathbf{q}$. In other words, a long-range magnetic wave order with a period of 4$b$ induces a strain wave with half of its period, 2$b$. The displacements of the atoms in CrOCl are resolved by neutron scattering\cite{zhang2014temperature}. Our density functional theory (DFT) calculations also confirm the displacements of atoms in the $\uparrow \uparrow \downarrow \downarrow$ ground state with a period of 2$b$, as shown in Fig. \ref{F1}b.\par

Based on the above phenomena, we can regard the $\uparrow \uparrow \downarrow \downarrow$ stripy AFM phase as an antiferroelectric order consisting of the electric dipoles of the distorted Cl-Cr-O chains away from the inversion center (dashed boxes in Fig. \ref{F1}b). Therefore, ordered effective dipoles are formed as indicated by the open arrows in Fig. \ref{F1}b, while the electric dipole of each atom is provided in Supplementary Information Fig. S1. Both the contributions of ions and electrons to polarization are taken into account by the Born effective charge method (see Methods). The relative permittivity $\epsilon_r$ of the bulk CrOCl crystal increases sharply at $\sim$14 K (the $\epsilon_r-T$ curve in Fig. \ref{F1}c) with the decreasing temperature, in perfect agreement with the Néel temperature where the net moment drops (the $M-T$ curve in Fig. \ref{F1}c), evidencing the structural phase transition to be a product of ME coupling. More importantly, the layer thickness experiences a sudden expansion at the Néel temperature\cite{reuvekamp2014investigation}, which tends to result in a decrease in the measured capacitance as $C=\frac{\epsilon_r\epsilon_0 S}{d}$ (see Methods). As a result, the increase in $\epsilon_r$ should be attributed to the appearance of additional antiferroelectric dipoles. According to the previous report\cite{gu2022magnetic}, CrOCl transforms into the collinear $\uparrow \uparrow \uparrow \downarrow \downarrow$ phase under an external out-of-plane field of $\sim$4.5 T (the $M-H$ curve in Fig. \ref{F1}d). Our calculations show that in the $\uparrow \uparrow \uparrow \downarrow \downarrow$ phase, the crystal relaxes to an orthorhombic structure (see Supplementary Information Table S1), in which the atomic displacements and the electric dipoles show a period of 5$b$. Under an external electric field of 0.07 V/{\AA} (close to the experimental values of the tunneling measurements below), the $\uparrow \uparrow \downarrow \downarrow$ phase produces a net polarization of $\rm{116.5\ \upmu C/m^2}$, which is eight times larger than that of the $\uparrow \uparrow \uparrow \downarrow \downarrow$ phase ($\rm{14.42\ \upmu C/m^2}$). This suggests that the electric dipoles of the $\uparrow \uparrow \downarrow \downarrow$ phase can be tilted by the external field, while the $\uparrow \uparrow \uparrow \downarrow \downarrow$ phase is rather structural stable and can hardly be adjusted by the external field. Correspondingly, the measured $\epsilon_r$ undergoes a decrease upon the magnetic phase transition to the $\uparrow \uparrow \uparrow \downarrow \downarrow$ phase (Fig. \ref{F1}d). The magnetic order-induced electric dipoles were further verified by the pyroelectric measurements, where the small net polarization could arise from uncompensated antiferroelectric domains (see Supplementary Information Fig. S2).\par

Now that we have demonstrated the existence of spin-induced electric dipoles in bulk CrOCl, from a macroscopic perspective, we would expect the interactions between the electric and magnetic degrees of freedom and consequently the manipulation of the spin order by the electric field, or vice versa. To achieve a large electric field in ultra-compact nanodevices, we fabricated vertical tunneling devices based on single- to few-layer CrOCl flakes (Fig. \ref{F2}a). Cross-structured few-layer graphite stripes are used to contact the CrOCl tunneling layer, and the whole device is encapsulated by hexagonal boron nitride (\textit{h}-BN). The electric field in CrOCl can reach $\sim$0.1 V/{\AA}, where the spin-induced electric dipoles are expected to be polarized and have a substantial impact on the tunneling current.\par

The existence of the antiferroelectric dipoles can be verified by the hysteresis of the $I-V$ curves. As shown in Fig. \ref{F2}b, the tunneling currents in the sweep-up and the sweep-down processes deviate significantly from each other under the same applied voltage at a temperature of 2 K, manifesting as an obvious clockwise hysteresis. In sharp contrast, the $I-V$ curve at 20 K (above the Néel temperature) shows no sign of hysteresis. The electric hysteresis can be further viewed in the 2D map of temperature and bias voltage (Fig. \ref{F2}c), where the current polarization $\rho$ is defined as $\rho=(R_{\rm{up}}-R_{\rm{down}})/(R_{\rm{up}}+R_{\rm{down}})$. Apart from the hysteresis loop induced by extrinsic charge traps at low voltages\cite{kaushik2017reversible}, $\rho$ exhibits non-zero values below $\sim$12 K, perfectly coincides with the Néel temperature of exfoliated CrOCl\cite{gu2022magnetic}. Clockwise hysteresis results in a negative $\rho$ under positive bias and a positive $\rho$ under negative bias. The absolute value of $\rho$ decreases with increasing temperature and eventually becomes zero at the Néel temperature (Fig. \ref{F2}c and Fig. \ref{F2}d), implying the disappearance of the interplay between charge and spin order. From the $R-T$ curves of different sweeping processes at 5 V (Fig. \ref{F2}e), it can be seen that in the sweep-up process, the resistance drops significantly (entering a lower resistance state) below the Néel temperature, while in the sweep-down process, the resistance increases (entering a higher resistance state). The above phenomena indicate that the configuration of the dipole order dominates the measured tunneling resistance below the Néel temperature. The resistance at $-5$ V exhibits completely antisymmetric behavior, as shown in Supplementary Information Fig. S3.\par

As mentioned earlier, in the absence of adjustable electric dipoles, the $\uparrow \uparrow \uparrow \downarrow \downarrow$ phase under an external magnetic field is rather structurally stable. The resistances under different electrical sweeping processes show completely different magnetoresistance behaviors at the $\uparrow \uparrow \downarrow \downarrow$ phase (see details in Supplementary Information Fig. S6 and S7) and converge at the transition field (Fig. \ref{F2}f), which is in perfect agreement with our theory. The hysteresis of the $I-V$ curve in the $\uparrow \uparrow \downarrow \downarrow$ state can be observed in all CrOCl devices of varying thicknesses down to monolayer. Suspending the sweeping process at a specific voltage, the resistance value relaxes slowly with a time constant of several hours (see Supplementary Information Fig. S5), which is highly reproducible as shown in a single-layer CrOCl device (Fig. \ref{F2}g). Such splitting of resistances between different electric sweep directions is direct evidence of ME coupling in 2D CrOCl, as the emergence of antiferromagnetic order leads to additional resistance states. The ramping electric field acts as an electric excitation to tilt the electric dipoles, while the interactions between spin and charge tend to relax the system to the antiferroelectric ground state. Hysteresis and relaxation behavior are discussed in detail in Supplementary Information Fig. S4 and S5. In addition to the $I-V$ hysteresis, the structural modulation of the $\uparrow \uparrow \downarrow \downarrow$ phase under the electric field also leads to different magnetoresistance behavior with varying bias voltage, as shown in Supplementary Information Fig. S8 and S9. \par

Switching the tunneling resistance between different metastable states resulting from ME coupling by sweeping the electric field provides us with design principles for magnetoelectronic devices. We have demonstrated that the magnetic transitions in CrOCl under external fields are first-order transitions with large hysteresis loops\cite{gu2022magnetic}. Therefore, we can expect that the electric field can adjust the tunneling resistance to any meta-stable states inside the magnetic hysteresis loop by tilting the electric dipoles. From the $I-B$ curve of the CrOCl tunneling device at 5.5 V (Fig. \ref{F3}a), the hysteresis closes at $B_s$=5.63 T, where CrOCl is believed to transform to the $\uparrow \uparrow \uparrow \downarrow \downarrow$ state. Subtracting the current of the $B$-up curve from that of the $B$-down curve, we obtain the maximum hysteresis occurs at $B_0$=3.33 T (Fig. \ref{F3}a). In other words, the spin configuration of CrOCl at 3.33 T is close to the $\uparrow \uparrow \downarrow \downarrow$ AFM ground state (low current) in the $B$-up cycle, but remains in the $\uparrow \uparrow \uparrow \downarrow \downarrow$ state (high current) in the $B$-down cycle. The maximum value of $\delta I$ is $\sim$2 nA, 10\% of the tunneling current. Similar to the operation in Fig. \ref{F2}, we apply electric excitations by sweeping the bias voltage up and back to 5.5 V. After the electric excitation, the tunneling current changes to an intermediate state following a slight relaxation and stabilizes at the state for several hours without any sign of change. The current value of the intermediate state is determined by the electric excitation, specifically, by the sweeping rate and the peak voltage. By applying different electric excitations by design, we can switch the tunneling current to an arbitrary value between the highest ($\uparrow \uparrow \uparrow \downarrow \downarrow$) and lowest ($\uparrow \uparrow \downarrow \downarrow$) current values, as demonstrated in Fig. \ref{F3}b. The noise fluctuation of each state is less than 0.02 nA, which means that a single tunneling junction can store at least a decimal number if the difference between distinguishable adjacent states is set to be an order of magnitude larger than the noise fluctuation. After sweeping the magnetic field across $B_s$ and back to $B_0$, the junction can be reset to the $\uparrow \uparrow \uparrow \downarrow \downarrow$ high current state, which means the erasing of the stored information.\par

Likewise, to expand the capabilities of the device, the magnetic field can also serve as another degree of freedom to tune the tunneling resistance. Sweeping the magnetic field to a value below $B_s$ and back to $B_0$ also adjusts the tunneling current to an intermediate state, as shown in Fig. \ref{F3}c. Accordingly, we can also adjust the tunneling current to different values with varying magnetic excitations, and set the current back to the $\uparrow \uparrow \downarrow \downarrow$ low current state with a large electric excitation. Repeating the previous electric excitations following different magnetic excitations, we obtain a new set of current values. In this way, we constructed a 5$\times$5 2D list of different resistance states by electric and magnetic excitations (Fig. \ref{F3}d). The opposite dependence of the current on electric and magnetic excitations strongly evidences the ME coupling in this system, as in the classical electrodynamics picture, the magnetic field flips the electron spins and the electric field reverses the electric dipoles. The combination of spin and charge order determines the tunneling resistance of the device, which can be manipulated both by electric and magnetic excitations. Furthermore, it is worth mentioning that similar multi-sate data storage operation principles are also applicable at zero-field and to single-layer CrOCl devices (see Supplementary Information Fig. S10 and S11). We finally verified the repeatability of device operation as shown in Fig. \ref{F3}e-f. The device is first set to the $\uparrow \uparrow \downarrow \downarrow$ AFM ground state ($\sim$ 19.4 nA) by a sufficiently large electric excitation of $\delta V= $2 V, and then reset to the $\uparrow \uparrow \uparrow \downarrow \downarrow$ state ($\sim$21.4 nA) by the magnetic stimulation of $\delta B= $1.67 T. The writing and erasing operations were repeated seven times, and the high and low resistance states showed perfect stability.\par

In summary, we have demonstrated the ME effect in the 2D stripy antiferromagnet CrOCl. By sweeping the bias voltage, we achieved electrically driven magnetic transitions of the AFM states. In particular, the tunneling resistance can be set to arbitrary values \textit{via} electric and magnetic excitations, which has never been reported in other 2D systems. Thanks to the stability of CrOCl, direct coupling between electric and magnetic degrees of freedom can be observed in all devices of varying thicknesses down to monolayer. Furthermore, the special ME coupling term gives rise to successive metastable states that are rather stable and hardly degenerate over time. The multi-state data storage realized in CrOCl may serve as a new paradigm with the potential to impact information technology, such as analog data storage and computation in an array of tunneling devices, stepping beyond Von Neumann architecture and enabling neuromorphic computing with low power consumption. Furthermore, the ME coupling in CrOCl may give rise to more unexplored fantastic properties, highlighting the characteristics of 2D antiferromagnetic materials and their promising potential in fundamental research and spintronic applications.\par

\vspace{-3mm}
\section*{Methods}
\vspace{-3mm}
\small{\textbf{Crystal synthesis and characterizations.}
The mixture of powdered $\rm{CrCl_3}$ and $\rm{Cr_2O_3}$ with a molar ratio of 1:1 and a total mass of 1.5 g was sealed in an evacuated quartz ampule. The ampule was then placed in a two-zone furnace where the source and sink temperatures for growth were set to 940 °C and 800 °C, respectively, and kept for two weeks. Subsequently, the furnace was slowly cooled to room temperature, and high-quality CrOCl crystals were obtained. Magnetization measurements were performed by standard modules of a Quantum Design PPMS.
}

\vspace{3mm}

\noindent \small{\textbf{Dielectric and pyroelectric measurements.} 
Both pyroelectric measurements and dielectric measurements were performed on a TeslatronPT System, Oxford Instruments. Silver epoxy was painted on opposite sides of the sample as electrodes. For pyroelectric measurements, an electrometer (Keithley 6517B) was used as a DC voltage source and ammeter. When the temperature stabilized at 60 K, an external electric field was applied to the sample. The sample was then cooled from 60 K to 2 K under different electric fields. After the temperature stabilized at 2 K, the electric field was removed and the temperature was increased to 60 K at a rate of 5 K/min. During the heating process, the change of pyroelectric current with temperature was collected. The polarization can be obtained by integrating the pyroelectric current with time. Dielectric measurements were made using a capacitive bridge meter (AH2700A). From 2 K to 50 K, the change in relative permittivity with temperature was measured by heating at the rate of 2 K/min. Magneto-dielectric effect at 2 K was measured at a magnetic field sweep rate of 15 Oe/s from 0 to 8 T. All test frequencies are 20 kHz.}

\vspace{3mm}

\noindent \small{\textbf{Device fabrication.}
Few-layer graphite, \textit{h}-BN (10-30 nm), and CrOCl flakes were obtained by the scotch tape exfoliation method under ambient conditions. The heterostructures were then assembled with a conventional tear-and-stack technique based on polypropylene carbonate(PPC)/polydimethylsiloxane(PDMS) polymer stacks placed on glass slides. Once encapsulated, the devices were annealed in a high vacuum with a mixed gas flow of $\rm{H_2}$ and Ar to remove residual PPC. Metal contact of Cr/Au (5/25 nm ) electrodes were then defined using electron-beam lithography, reactive ion etching (in plasma of the $\rm{CHF_3/O_2}$ mixture), electron beam evaporation, and lift-off processes.}

\vspace{3mm}

\noindent \small{\textbf{Electrical transport measurements.} 
Transport measurements were performed in a $\rm{Heliox^3}$ He insert system equipped with a 14 T superconducting magnet. The lowest temperature of the system is 1.6 K. To measure the $I-V$ characteristics of the tunnel barrier and the magnetoresistance, a Keithley 2636B source meter was used to apply a bias voltage and a standard two-probe module was used to measure the tunneling current. To obtain intrinsic signals and at the same time exclude the possibility of the Joule heating effect, the tunneling current is limited to $\sim$50 nA, so the total power in a junction of $\sim$1 ${\upmu}\rm{m^2}$ is merely $\sim$0.3 $\upmu$W.}


\vspace{3mm}

\noindent \small{\textbf{DFT calculations.} 
Our DFT calculations were performed using the generalized gradient approximation for the exchange-correlation potential, the projector augmented wave method\cite{kresse1999ultrasoft}, and a plane-wave basis set implemented in the Vienna ab-initio simulation package (VASP)\cite{kresse1996efficiency}. Dispersion correction was made at the van der Waals density functional (vdW-DF) level\cite{dion2004van}, with the optB86b functional for the exchange potential\cite{klimevs2011van}, and which was proved to be accurate in describing the structural properties of layered materials\cite{qiao2018few} and was adopted for structure related calculations. The shape and volume of each supercell and all atomic positions were fully relaxed until the residual force per atom was less than $1\times 10^{-3} $ eV\cdot {\AA}$^{-1}$ in our calculations. In VASP calculations, the kinetic energy cut-off for the plane-wave basis set was set to be 700 eV for geometric and electronic structure calculations. A $k$-mesh of $10\times14\times4$ was adopted to sample the first Brillouin zone of the conventional unit cell of CrOCl bulk. The $U$ and $J$ values of the on-site Coulomb interaction of the Cr $d$ orbitals are 3.0 eV and 1.0 eV, respectively, as revealed by a linear response method\cite{cococcioni2005linear} and comparison with the experimental results\cite{gu2022magnetic}. These values are comparable to those adopted in modeling CrSCl\cite{wang2019family} and $\rm{CrI_3}$\cite{jiang2019stacking}. The Born effective charges were calculated using density functional perturbation theory\cite{baroni2001phonons}. The dipole moment of each atom is calculated by $P_i=Z_i^*\cdot \mu_i$, where $P_i$ is the dipole moment of the ion $i$ in one unit cell, $Z_i$ is the Born effective charge tensors and $\mu_i$ is the atomic displacement. For calculations of single-layer CrOCl, a sufficiently large vacuum layer over 20 {\AA} along the out-of-plane direction was adopted to eliminate the interaction among monolayer. The out-of-plane electric polarization of single-layer CrOCl under the external electric field is well defined in terms of the classical model due to the presence of a vacuum layer, which is calculated by integrating electron density times $z$-coordinate over the supercell.}

\vspace{3mm}

\section*{Data Availability}
\normalsize{The data that support the findings of this study will be available at an open-access repository with a doi link, when accepted for publishing.}

\section*{Code Availability}
\normalsize{The codes that support the findings of this study will be available at an open-access repository with a doi link, when accepted for publishing.}

\section*{Acknowledgements} 
\normalsize{This work was supported by the National Key R\&D Program of China (Grants No. 2018YFA0306900), and Beijing Natural Science Foundation (Grant No. JQ21018). Y.S. acknowledge support from the National Natural Science Foundation of China (Grant No. 51725104). W.J. acknowledge support from the Ministry of Science and Technology (MOST) of China (Grant No. 2018YFE0202700), and the National Natural Science Foundation of China (Grants No. 61761166009, No. 11974422 and No. 12104504). C.W. was supported by the China Postdoctoral Science Foundation (2021M693479). Calculations were performed at the Physics Lab of High-Performance Computing of Renmin University of China, Shanghai Supercomputer Center. K.W. and T.T. acknowledge support from the Elemental Strategy Initiative conducted by the MEXT, Japan (Grant NO. JPMXP0112101001), JSPS KAKENHI (Grant Nos. 19H05790, 20H00354 and 21H05233) and the A3 Foresight by JSPS.}



\section*{Authors contributions} 
\normalsize{Y.Y. and P.G. conceived the project. P.G. synthesized the CrOCl crystals and fabricated the devices. P.G. conducted the transport measurements with the help of Z.D., Q.W. and Z.H. C.W. conducted the DFT calculations under the supervision of W.J. D.S. conducted the dielectric and pyroelectric measurement under the supervision of Y.S. K.W. and T.T. grew the \textit{h}-BN bulk crystals. P.G. and Y.Y. drafted the manuscript. All authors discussed the results and contributed to the manuscript. }

\section*{Competing interests} 
\normalsize{The authors declare no competing interests.}

\bibliographystyle{naturemag}
