††thanks: These authors contributed equally.††thanks: These authors
contributed equally.

# Deterministic and Efficient Switching of Sliding Ferroelectrics

Shihan Deng  Key Laboratory of Computational Physical Sciences (Ministry of
Education), Institute of Computational Physical Sciences, State Key Laboratory
of Surface Physics, and Department of Physics, Fudan University, Shanghai
200433, China  Shanghai Qi Zhi Institute, Shanghai 200030, China  Hongyu Yu
Key Laboratory of Computational Physical Sciences (Ministry of Education),
Institute of Computational Physical Sciences, State Key Laboratory of Surface
Physics, and Department of Physics, Fudan University, Shanghai 200433, China
Shanghai Qi Zhi Institute, Shanghai 200030, China  Junyi Ji  Key Laboratory of
Computational Physical Sciences (Ministry of Education), Institute of
Computational Physical Sciences, State Key Laboratory of Surface Physics, and
Department of Physics, Fudan University, Shanghai 200433, China  Shanghai Qi
Zhi Institute, Shanghai 200030, China  Changsong Xu<EMAIL_ADDRESS>Key
Laboratory of Computational Physical Sciences (Ministry of Education),
Institute of Computational Physical Sciences, State Key Laboratory of Surface
Physics, and Department of Physics, Fudan University, Shanghai 200433, China
Shanghai Qi Zhi Institute, Shanghai 200030, China  Hongjun Xiang
<EMAIL_ADDRESS>Key Laboratory of Computational Physical Sciences
(Ministry of Education), Institute of Computational Physical Sciences, State
Key Laboratory of Surface Physics, and Department of Physics, Fudan
University, Shanghai 200433, China  Shanghai Qi Zhi Institute, Shanghai
200030, China

###### Abstract

Recent studies highlight the scientific importance and broad application
prospects of two-dimensional (2D) sliding ferroelectrics, which prevalently
exhibit vertical polarization with suitable stackings. It is crucial to
understand the mechanisms of sliding ferroelectricity and to deterministically
and efficiently switch the polarization with optimized electric fields. Here,
applying our newly developed DREAM-Allegro multi-task equivariant neural
network, which simultaneously predicts interatomic potentials and Born
effective charges, we construct a comprehensive potential for boron nitride
($\mathrm{BN}$) bilayer. The molecular dynamics simulations reveal a
remarkably high Curie temperature of up to 1500K, facilitated by robust
intralayer chemical bonds and delicate interlayer van der Waals(vdW)
interactions. More importantly, it is found that, compared to the out-of-plane
electric field, the inclined field not only leads to deterministic switching
of electric polarization, but also largely lower the critical strength of
field, due to the presence of the in-plane polarization in the transition
state. This strategy of an inclined field is demonstrated to be universal for
other sliding ferroelectric systems with monolayer structures belonging to the
symmetry group $p\bar{6}m2$, such as transition metal dichalcogenides (TMDs).

Achieving ferroelectric devices with small thickness and perpendicular
polarization is a critical step toward realizing low-energy, nonvolatile, and
high-density ferroelectric memory [1]. However, it is well-known that the
depolarization field becomes significant in thin film ferroelectrics.
Additionally, “dead layers” inevitably form due to factors such as interface
effects, defects, and impurities, all of which lead to performance degradation
in the application of ultrathin films [2, 3, 4, 5].

The recently proposed concept of sliding ferroelectricity in 2D vdW materials
naturally overcomes these obstacles [6]. The prototype of sliding
ferroelectrics is $\mathrm{BN}$ bilayers [see Fig. 1(a)], which display energy
ground state as the AB (BA) configuration [7, 2, 1, 8, 9, 10] with $C_{3v}$
symmetry. Within the AB configuration, the top $\mathrm{B}$ atom aligns with
the bottom $\mathrm{N}$ atom, and the top $\mathrm{N}$ atom and the bottom
$\mathrm{B}$ atom sit at the hollow sites of the honeycomb lattice. The AB(BA)
stacking mode disrupts mirror symmetry and leads to an out-of-plane
polarization of 2 pC/m in the $-z$($z$) direction [6, 11]. This vertical
polarization arises from charge redistribution driven by interlayer coupling
[2, 12, 4, 13]. The BA state can transition to the AB state via relative
sliding between the two layers by a $\mathrm{B}$-$\mathrm{N}$ bond length in
any of the three rotationally symmetric directions [11, 13, 14]. Conversely,
sliding in the opposite direction by one bond length leads to the mirror-
symmetric ($D_{3h}$) AA state, representing the highest energy state. Sliding
ferroelectrics, especially the proposal of bilayer stacking ferroelectricity
(BSF) theory [15], broaden the spectrum of candidate materials for 2D
ferroelectrics, as bilayers made of nonpolar monolayers can exhibit
ferroelectricity through specific stackings. Beyond $\mathrm{BN}$ bilayer,
sliding ferroelectrics has been experimentally confirmed in various vdW
systems, including semimetals [16, 17, 18, 19, 20], insulators [21],
semiconductors [22, 23, 24, 25, 26], and organic crystals [27], all of which
are robust at room temperature. Notably, recent experiments highlighting high
endurance and fatigue resistance in sliding ferroelectrics underscore the
significant potential of these materials for practical applications [28, 29].

The study of sliding ferroelectricity is rapidly developing, with both
fundamental and practical challenges yet to be addressed. Although the
switching energy barriers are typically very low ($\sim$meV per unit cell [6,
30, 31]), sliding ferroelectrics exhibit high Curie temperatures ($T_{C}$), as
most of them are operated at room-temperature. Thermodynamic modeling using a
mean-field approximation suggests $T_{C}=1.58\times 10^{4}$ K for
$\mathrm{BN}$ bilayer [14]. Such prediction reflects the good stability of
$\mathrm{BN}$ bilayers. Moreover, the polarization of sliding ferroelectrics
is typically small, leading to a large coercive field for polarization
switching, which indicates that a wiser and more efficient switching strategy
is highly desired. Symmetry changes during slidings indicating that the Born
effective charges (BECs) can change dramatically [32], the treatment of fixed
BEC may be insufficient [33] [see Fig.S1 in the supplemental material (SM)
[34] for further discussion]. Ab initio molecular dynamics (AIMD) can be used
to study sliding ferroelectricity, but it cannot handle large systems and is
computationally consuming [51]. Hence, an accurate and efficient description
of both intralayer and interlayer couplings is required.

Figure 1:  (a) Schemes of the AA, BA, AB, and SP stackings of $\mathrm{BN}$
bilayers. Blue arrows define the lattice vectors $\bm{a}$, $\bm{b}$ and
$\bm{c}$ of unit cell. Purple arrows indicate the sliding vector, $\bm{t}$,
which is defined as the fractional displacement of the top layers with respect
to bottom layer in the AA stacking. Energy $\varepsilon_{0}$ (b) and
polarization $P$ (c) along the designed sliding pathway, as calculated by the
DFT and Berry phase approach [52, 53].

In this Letter, applying our newly developed DREAM method (i.e., a generalized
dielectric response equivariant atomistic multi-task framework based on
Allegro [45]), we construct a neural network model for $\mathrm{BN}$ bilayer.
Such model is capable of predicting both the interatomic potential and the BEC
tensors, which enable prediction of a reasonable $T_{C}$ and accurate
responses to applied electric fields. Moreover, it is found that an inclined
electric field, which breaks the three-fold rotational symmetry of
$\mathrm{BN}$ bilayer, can not only deterministically control directions of
slidings, but also substantially reduce the total coercive field. These
findings are general and can be applied to other similar sliding
ferroelectrics.

blueFerroelectric transition of BN bilayer. Firstly, the basic properties of
$\mathrm{BN}$ bilayers are examined by density functional theory (DFT)
calculations. As illustrated in Fig. 1(b)(c), the out-of-plane polarization
($P_{\perp}$) of the AB (BA) state yields $\mp$2.078 pC/m, and the energy
barrier for polarization switching is determined to be 7.44 meV per unit cell.
Such results are well consistent with those of previous studies [6, 11].
Notably, the transition state or intermediate saddle-point (SP) state with
$C_{2v}$ symmetry exhibits a strictly in-plane polarization ($P_{\parallel}$)
of 2.091 pC/m.

To train the potential with the original Allegro method [45], we run on-the-
fly machine learning force fields(MLFFs) [40, 41, 42] and generate 1659
structures with $5\times 5\times 1$ supercell, across a wide temperature range
below 3000 K. The structure, energy, force, and stress tensor are included in
training. The obtained model can accurately predict the energy difference
among structural configurations, with a small mean absolute errors (MAEs) of
0.053 meV/atom [see SM [34] for details].

Figure 2:  The change of $\mathrm{Im}(\psi(\bm{t}))$ as a function of
temperature. Red dots and blue triangles represent the heating and cooling
processes, respectively, which are marked by red and blue arrows. Green arrows
denote the reciprocal lattice vectors $\bm{a}^{*}$ and $\bm{b}^{*}$, and the
black arrow indicates that $\bm{G}=\bm{a}^{*}-\bm{b}^{*}$ is selected. Purple
arrows visualize the sliding vector $\bm{t}$ of the BA and AB state, which are
$\bm{t}_{\mathrm{BA}}=1/3(\bm{b}-\bm{a})$ and
$\bm{t}_{\mathrm{AB}}=2/3(\bm{b}-\bm{a})$ respectively. The dashed gray line
signifies the critical temperature.

In the case of sliding ferroelectricity, $T_{C}$ can be characterized by the
sliding vector $\bm{t}$ [14]:

$\lim\limits_{T\to T_{C}^{-}}\frac{\partial\langle\bm{t}\rangle}{\partial
T}\to\infty,$ (1)

where $\langle~{}\rangle$ signifies the ensemble average at a certain
temperature. Vector $\bm{t}$ is the averaged displacements of the top layer
with respect to the bottom layer in the AA stacking [see Fig. 1(a)]. To avoid
the impact of in-plane periodic boundary conditions, we use the reciprocal
lattice vector ($\bm{G}$) to define a new order parameter ($\psi$), which
reads

$\psi(\bm{t})=e^{i\bm{G}\cdot\bm{t}}.$ (2)

where $\psi$ ranges between $[-1,1]$ and remains unchanged under the
substitution $\bm{t}\to\bm{t}+\bm{R}$, where $\bm{R}$ represents the in-plane
lattice vector.

We now focus on the ferroelectric-paraelectric transition of $\mathrm{BN}$
bilayer. The MD simulations are performed using our potential model, where 900
atoms is incrementally heated and cooled between 1 K and 2000 K, with an
interval of 50 K. At each temperature, the quantity of $\psi$ is computed for
the equilibrium states over 50 ps. The changes in the imaginary part of $\psi$
as a function of temperature are presented in Fig. 2 for
$\bm{G}=\bm{a}^{*}-\bm{b}^{*}$. Starting in the BA state
[$\mathrm{Im}(\psi_{\mathrm{BA}})=\sqrt{3}/2$], an increase in temperature
progressively triggers system sliding. Conversely, during the cooling process,
system sliding gradually ceases and eventually randomly stabilizes in the AB
state [$\mathrm{Im}(\psi_{\mathrm{AB}})=-\sqrt{3}/2$]. Accoriding to Fig. 2,
the critical temperature is approximated as 1500 K. Intriguingly, even for
temperatures over $T_{C}$, the system conducts only few complete slides within
50 ps, oscillating around the BA or AB state most of the time. This phenomenon
causes $\psi$ to oscillate above and below zero, contrasting with traditional
order parameters that promptly vanish above $T_{C}$. The presently predicted
$T_{C}=1500$ K is much more reasonable than the $1.58\times 10^{4}$ K from
mean-field approximation [14].

blueSwitching of sliding ferroelectric $\mathrm{BN}$ bilayer. To simulate the
$\mathrm{BN}$ bilayer under an external electric field, we integrate Allegro
within the new DREAM framework [46, 34]. This integration enables the network
to learn both scalars and tensors by leveraging geometric tensors [54, 55].
Consequently, the network is capable of simultaneously predicting interatomic
potentials and BECs. The MAEs of the final model, which incorporates BEC
information of 3484 structures at temperatures below 500K, yields 0.030
meV/atom in energy and 0.008 $|$e$|$ in BEC, demonstrating good accuracy of
the model.

Figure 3:  (a) Energy $\varepsilon_{0}$ and polarization $P$ distribution
calculated by the DFT and Berry phase approach. Positions within the hexagon
correspond to different sliding vectors, with the center representing the BA
state. The shades of grey in the background indicate energy levels. The
direction of colored arrows shows in-plane polarizations $P_{\parallel}$,
while their colors represent out-of-plane polarization $P_{\perp}$. The
sliding directions $\theta$ are indicated by large hollow arrows. (b) The
energy surface with only $E_{\perp}$ along $-z$ (the top panel) and with both
$E_{\perp}$ along $-z$ and $E_{\parallel}$ at $60^{\circ}$ (the bottom panel),
calculated according to Eq. 3 with $E_{\perp}$ and $E_{\parallel}$ equal to
1V/Å. White arrows denote the sliding pathways with the lowest barrier. (c)
Energy along the sliding pathway under varying $E_{\parallel}$ and
$E_{\perp}$, also calculated using Eq. 3. The $E_{\parallel}$ increases from
left to right panels, while $E_{\perp}$ are gradually enhanced from top to
bottom panels.

We now examine the responses of the $\mathrm{BN}$ bilayer to perpendicular
field $E_{\perp}$. We start with the ferroelectric BA state and incrementally
increase $E_{\perp}$ along the $-z$ direction over time at a rate of 0.02
V/(Å$\cdot$ps), at temperatures ranging from 100 K to 500 K. The sliding
process occurs rapidly, typically within a few picoseconds, and thermal
vibrations add complexity to identifying the exact onset of sliding.
Therefore, we define the coercive field ($E_{\perp,c}$) as the field at which
sliding reaches the SP state. Our results indicate that $E_{\perp,c}$
decreases with increasing temperature, starting at approximately 1.99 V/Å at
100 K and saturating at around 1.39 V/Å for over 200 K. Moreover, since
$E_{\perp}$ does not break the three-fold rotational symmetry of the
$\mathrm{BN}$ bilayer, sliding occurs randomly along the three rotationally
symmetric directions in each simulation [see the top panel of Fig. 3(b)].

Given that some sliding states of $\mathrm{BN}$ bilayer exhibit
$P_{\parallel}$, we investigate the effect of in-plane electric field
($E_{\parallel}$) on ferroelectric switching. As revealed by Fig. 3(a), the
direction of $P_{\parallel}$ is the same as the corresponding sliding
direction. These directions can be characterized by the angle $\theta$, which
is measured through a counterclockwise rotation from the $y$-axis. This
finding implies that $E_{\parallel}$ oriented along one of these three
directions can break the three-fold rotational symmetry and induce
deterministic switching. For instance, $E_{\parallel}$ oriented at
$60^{\circ}$ will lead to sliding being only along $60^{\circ}$ [see the lower
panel of Fig. 3(b)]. Actually, applying $E_{\parallel}$ within the range of
$(0^{\circ},120^{\circ})$ can all leads to sliding along $60^{\circ}$ [refer
to Fig.S4]. In addition, Fig. 3(a) shows significant $P_{\parallel}$ at the
bridging SP state, which indicates that suitable $E_{\parallel}$ can
substantially reduces the energy barrier. Such conjecture can be verified by
energy distribution shown in the first two columns of Fig. 3(c)], where the
energy barrier for sliding from BA to AB state is obviously reduced in the
presence of $E_{\parallel}$.

Figure 4:  (a) $E_{\perp,c}$ as a function of temperature. The orange and
purple lines represent the scenarios with $E_{\parallel}=0$ and
$E_{\parallel}=0.2$V/Å, respectively, at $\theta=60^{\circ}$. The error bars
indicate the range of results from parallel simulations. (b) The relationship
between $E_{\perp,c}$ and $E_{\parallel}$, as obtained by analytical
calculations (purple line), numerical calculations (green line), and MD
simulations at 0.1K (yellow line). The pentagram marks on the curves indicate
the the minimum $E_{t,c}$. The insets depict the energy along the sliding
pathway with the minimum $E_{t,c}$ obtained by analytical (the top panel) and
numerical (the bottom panel) calculations, respectively. (c) The magnitude
(blue line) and dip angle (red line) of the minimum $E_{t,c}$ across different
temperatures. The dip angle is defined as the angle between $E_{t,c}$ and the
horizontal plane.

Moreover, the decreasing in the energy barriers also indicate that the
presence of $E_{\parallel}$ may reduce $E_{\perp,c}$. To examine such
conjecture, we fit the energy $\varepsilon_{0}$, in-plane dipole
$p_{\parallel}$, and out-of-plane dipole $p_{\perp}$ as a function of the
amount of sliding $t$ with respect to the SP state, along the BA to AB sliding
pathway. The fitted models yield
$\varepsilon_{0}(t)=e_{4}t^{4}-2e_{4}t_{\mathrm{BA}}^{2}t^{2}+e_{0}$,
$p_{\parallel}(t)=p_{4}^{\parallel}t^{4}+p_{2}^{\parallel}t^{2}-p_{4}^{\parallel}t_{\mathrm{BA}}^{4}-p_{2}^{\parallel}t_{\mathrm{BA}}^{2}$,
and $p_{\perp}(t)=p_{3}^{\perp}t^{3}-3p_{3}^{\perp}t_{\mathrm{BA}}^{2}t$,
where $t_{\mathrm{BA}}=-1/6|\bm{b}-\bm{a}|=-0.725$ Å is a constant. The
coefficients are determined to be $e_{4}=27.017$ meV/Å4, $e_{0}=7.41$ meV,
$p_{4}^{\parallel}=5.82\times 10^{-3}$ e/Å3, $p_{2}^{\parallel}=-1.66\times
10^{-2}$ e/Å, and $p_{3}^{\perp}=9.13\times 10^{-3}$ e/Å2 [refer to Fig.S3].
The total potential energy after applying $E_{\parallel}$ along the BA to AB
sliding pathway and $E_{\perp}$ can then be expressed as

$\varepsilon(t)=\varepsilon_{0}(t)-E_{\parallel}\cdot
p_{\parallel}(t)-E_{\perp}\cdot p_{\perp}(t).$ (3)

If we define the critical field as the minimum field at which the local energy
minimum of $\varepsilon(t)$ disappears [i.e.,
$\mathrm{d}\varepsilon(t)/\mathrm{d}t\leq 0$, see the top inset of Fig. 4(b)]
along the BA to SP sliding pathway. Then, $E_{\perp,c}$ along $-z$ can be
obtained analytically for a given $E_{\parallel}$ by solving Eq. 3 [see the
purple line in Fig. 4(b) and SM [34]]. When $E_{\parallel}$ is small,
$E_{\perp,c}$ can be expressed in terms of $\sqrt{E}_{\parallel}$ as

$\displaystyle E_{\perp,c}=$
$\displaystyle-\frac{4e_{4}t_{\mathrm{BA}}}{3p_{3}^{\perp}}\left(1-\sqrt{-\frac{p_{2}^{\parallel}+2p_{4}^{\parallel}t_{\mathrm{BA}}^{2}}{e_{4}t_{\mathrm{BA}}^{2}}}\sqrt{E_{\parallel}}\right)$
(4) $\displaystyle=$ $\displaystyle
E^{0}_{\perp,c}(1-\beta\sqrt{E_{\parallel}}),$

where $\beta>0$. Clearly, Eq. 4 shows that the presence of $E_{\parallel}$
will reduce the required $E_{\perp,c}$. In the absence of $E_{\parallel}$,
$E^{0}_{\perp,c}=-4e_{4}t_{\mathrm{BA}}/3p_{3}^{\perp}=$ 2.86 V/Å is the
vertical coercive field. More importantly, it is further found that the
presence of $E_{\parallel}$ can lower the total coercive field $E_{t,c}$. When
$E_{\parallel}\to 0$, $E_{t,c}$ can be expressed as

$\displaystyle E_{t,c}$
$\displaystyle=\sqrt{(E_{\perp,c})^{2}+(E_{\parallel})^{2}}$ (5)
$\displaystyle\to
E^{0}_{\perp,c}\left(1-\beta\sqrt{E_{\parallel}}\right)<E^{0}_{\perp,c},$

which clearly demonstrates that $E_{t,c}$ is smaller than $E^{0}_{\perp,c}$
with finite $E_{\parallel}$. Alternatively, if we define the critical field as
the minimum field such that $\varepsilon(t)$ does not exceed that at BA [i.e.,
$\varepsilon(t)\leq\varepsilon(t_{\mathrm{B}A})$, see the bottom inset of Fig.
4(b)], $E_{\perp,c}$ can also be solved numerically for a given
$E_{\parallel}$, as shown by the green line in Fig. 4(b). In both analytical
and numerical cases, the minimum $E_{t,c}$ is reduced by more than 64% and
70%, respectively, comparing to the $E^{0}_{\perp,c}$.

Then, we conduct MD simulations to verify the prediction that finite
$E_{\parallel}$ can lower $E_{\perp,c}$ and $E_{t,c}$. By applying a constant
$E_{\parallel}$ of 0.2 V/Å along $\theta=60^{\circ}$ direction, $E_{\perp,c}$
is found to significantly reduce, comparing to $E_{\parallel}=0$ [see Fig.
4(a)]. As in the previous analysis, this adjustment ensured that sliding
occurred exclusively at $\theta=60^{\circ}$. However, insufficient
$E_{\parallel}$ magnitude may not adequately counteract random thermal
fluctuations, potentially resulting in alternative sliding directions.

To determine the configurations of $E_{t,c}$ necessary to initiate sliding, we
apply an inclined electric field to the $\mathrm{BN}$ bilayer at different
temperatures. The field has varying in-plane and vertical components, adjusted
in increments of 0.05 V/Å, to initiate sliding within 30 ps [refer to Fig.S5].
At 0.1 K, the configuration of $E_{t,c}$ is shown by the yellow line in Fig.
4(b), exhibiting a trend consistent with previous theoretical analyses. As
shown in Fig. 4(c), the minimum $E_{t,c}$ decreases by over 70% as the
temperature decreases, compared to the vertical $E_{t,c}$. It is noteworthy
that experimentally measured coercive fields are often significantly lower
than theoretically predicted values. This discrepancy, potentially related to
the Landauer paradox [56, 57, 58], may be due to the presence of
inhomogeneities in the experimental samples [59, 60]. On another aspect, dip
angle of the minimum $E_{t,c}$ at different temperatures almost keep constant
of $24^{\circ}$. Such constant behavior will benefit to experimental setups
and also to device design.

Discussion. The decrease in $E_{t,c}$ with applying an inclined electric field
relies on the presence of in-plane polarization along the sliding pathway.
According to symmetry analysis [15], bilayers, which are made of monolayers
with the same layer group of $p\bar{6}m2$ as $\mathrm{BN}$, all display the
required in-plane polarization. Such discovery indicates that the strategy of
inclined field can be generalized to different sliding systems, such as widely
studied transition metal dichalcogenides (TMDs) $\mathrm{WSe_{2}}$,
$\mathrm{MoSe_{2}}$, $\mathrm{WS_{2}}$, and $\mathrm{MoS_{2}}$ [24, 25, 29],
as well as $\mathrm{InSe}$ [22, 23]. To confirm this prediction, we verify the
energy and polarization distribution of $\mathrm{MoS_{2}}$ [refer to Fig.S6],
which differs from $\mathrm{BN}$ only in the opposite direction of in-plane
polarization. Hence, the direction of $E_{\parallel}$ need to be reversed to
achieve efficient switching.

In conclusion, our newly developed machine-learning-assisted DREAM-Allegro
model effectively predicts the interatomic energy and BECs of parallel-stacked
$\mathrm{BN}$ bilayers. This enables accurate prediction on temperature and
field effects by performing MD simulations. Applying both simulations and
theoretical analysis, we find that in-plane field can lead to deterministic
switching of polarization and an inclined field can largely reduce the
critical field required for switching. Such strategy is further generated to
widely studied sliding ferroelectrics and thus paves the way for applications.

###### Acknowledgements.

We acknowledge financial support from the National Key R&D Program of China
(No. 2022YFA1402901), NSFC (No. 11825403, 11991061, 12188101, 12174060, and
12274082), the Guangdong Major Project of the Basic and Applied Basic Research
(Future functional materials under extreme conditions—2021B0301030005),
Shanghai Science and Technology Program (No. 23JC1400900), and Shanghai Pilot
Program for Basic Research—FuDan University 21TQ1400100 (23TQ017). C.X. also
acknowledges support from the Shanghai Science and Technology Committee (Grant
No. 23ZR1406600) and the Xiaomi Young Talents Program.

## References

  * Yasuda  _et al._ [2021] K. Yasuda, X. Wang, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Stacking-engineered ferroelectricity in bilayer boron nitride, Science 372, 1458 (2021), https://www.science.org/doi/pdf/10.1126/science.abd3230 .
  * Stern  _et al._ [2021] M. V. Stern, Y. Waschitz, W. Cao, I. Nevo, K. Watanabe, T. Taniguchi, E. Sela, M. Urbakh, O. Hod, and M. B. Shalom, Interfacial ferroelectricity by van der Waals sliding, Science 372, 1462 (2021), https://www.science.org/doi/pdf/10.1126/science.abe8177 .
  * Müller  _et al._ [2015] J. Müller, P. Polakowski, S. Mueller, and T. Mikolajick, Ferroelectric hafnium oxide based materials and devices: Assessment of current status and future prospects, ECS Journal of Solid State Science and Technology 4, N30 (2015).
  * Wang  _et al._ [2023] C. Wang, L. You, D. Cobden, and J. Wang, Towards two-dimensional van der Waals ferroelectrics, Nature Materials 22, 542 (2023).
  * Stengel and Spaldin [2006] M. Stengel and N. A. Spaldin, Origin of the dielectric dead layer in nanoscale capacitors, Nature 443, 679 (2006).
  * Li and Wu [2017] L. Li and M. Wu, Binary compound bilayer and multilayer with vertical polarizations: Two-dimensional ferroelectrics, multiferroics, and nanogenerators, ACS Nano 11, 6382 (2017).
  * Ribeiro and Peres [2011] R. M. Ribeiro and N. M. R. Peres, Stability of boron nitride bilayers: Ground-state energies, interlayer distances, and tight-binding description, Phys. Rev. B 83, 235312 (2011).
  * Wu and Li [2021] M. Wu and J. Li, Sliding ferroelectricity in 2d van der Waals materials: Related physics and future opportunities, Proceedings of the National Academy of Sciences 118, e2115703118 (2021), https://www.pnas.org/doi/pdf/10.1073/pnas.2115703118 .
  * Constantinescu  _et al._ [2013] G. Constantinescu, A. Kuc, and T. Heine, Stacking in bulk and bilayer hexagonal boron nitride, Phys. Rev. Lett. 111, 036104 (2013).
  * Tsymbal [2021] E. Y. Tsymbal, Two-dimensional ferroelectricity by design, Science 372, 1389 (2021), https://www.science.org/doi/pdf/10.1126/science.abi7296 .
  * Jiang  _et al._ [2022] W. Jiang, C. Liu, X. Ma, X. Yu, S. Hu, X. Li, L. A. Burton, Y. Liu, Y. Chen, P. Guo, X. Kong, L. Bellaiche, and W. Ren, Anomalous ferroelectricity and double-negative effects in bilayer hexagonal boron nitride, Phys. Rev. B 106, 054104 (2022).
  * Zheng  _et al._ [2020] Z. Zheng, Q. Ma, Z. Bi, S. de la Barrera, M.-H. Liu, N. Mao, Y. Zhang, N. Kiper, K. Watanabe, T. Taniguchi, J. Kong, W. A. Tisdale, R. Ashoori, N. Gedik, L. Fu, S.-Y. Xu, and P. Jarillo-Herrero, Unconventional ferroelectricity in moiré heterostructures, Nature 588, 71 (2020).
  * Liu  _et al._ [2023] M. Liu, H. Ji, Z. Fu, Y. Wang, J.-T. Sun, and H.-J. Gao, Orbital distortion and electric field control of sliding ferroelectricity in a boron nitride bilayer, Journal of Physics: Condensed Matter 35, 235001 (2023).
  * Tang and Bauer [2023] P. Tang and G. E. W. Bauer, Sliding phase transition in ferroelectric van der Waals bilayers, Phys. Rev. Lett. 130, 176801 (2023).
  * Ji  _et al._ [2023] J. Ji, G. Yu, C. Xu, and H. J. Xiang, General theory for bilayer stacking ferroelectricity, Phys. Rev. Lett. 130, 146801 (2023).
  * Fei  _et al._ [2018] Z. Fei, W. Zhao, T. A. Palomaki, B. Sun, M. K. Miller, Z. Zhao, J. Yan, X. Xu, and D. H. Cobden, Ferroelectric switching of a two-dimensional metal, Nature 560, 336 (2018).
  * Sharma  _et al._ [2019] P. Sharma, F.-X. Xiang, D.-F. Shao, D. Zhang, E. Y. Tsymbal, A. R. Hamilton, and J. Seidel, A room-temperature ferroelectric semimetal, Science advances 5, eaax5080 (2019).
  * Xiao  _et al._ [2020] J. Xiao, Y. Wang, H. Wang, C. D. Pemmaraju, S. Wang, P. Muscher, E. J. Sie, C. M. Nyby, T. P. Devereaux, X. Qian, X. Zhang, and A. M. Lindenberg, Berry curvature memory through electrically driven stacking transitions, Nature Physics 16, 1028 (2020).
  * de la Barrera  _et al._ [2021] S. C. de la Barrera, Q. Cao, Y. Gao, Y. Gao, V. S. Bheemarasetty, J. Yan, D. G. Mandrus, W. Zhu, D. Xiao, and B. M. Hunt, Direct measurement of ferroelectric polarization in a tunable semimetal, Nature Communications 12, 5298 (2021).
  * Jindal  _et al._ [2023] A. Jindal, A. Saha, Z. Li, T. Taniguchi, K. Watanabe, J. C. Hone, T. Birol, R. M. Fernandes, C. R. Dean, A. N. Pasupathy, and D. A. Rhodes, Coupled ferroelectricity and superconductivity in bilayer Td-MoTe2, Nature 613, 48 (2023).
  * Lv  _et al._ [2022] M. Lv, X. Sun, Y. Chen, T. Taniguchi, K. Watanabe, M. Wu, J. Wang, and J. Xue, Spatially resolved polarization manipulation of ferroelectricity in twisted hBN, Advanced Materials 34, 2203990 (2022), https://onlinelibrary.wiley.com/doi/pdf/10.1002/adma.202203990 .
  * Hu  _et al._ [2019] H. Hu, Y. Sun, M. Chai, D. Xie, J. Ma, and H. Zhu, Room-temperature out-of-plane and in-plane ferroelectricity of two-dimensional $\beta$-InSe nanoflakes, Applied Physics Letters 114, 252903 (2019), https://pubs.aip.org/aip/apl/article-pdf/doi/10.1063/1.5097842/14525668/252903_1_online.pdf .
  * Sui  _et al._ [2023] F. Sui, M. Jin, Y. Zhang, R. Qi, Y.-N. Wu, R. Huang, F. Yue, and J. Chu, Sliding ferroelectricity in van der Waals layered $\gamma$-InSe semiconductor, Nature Communications 14, 36 (2023).
  * Wang  _et al._ [2022] X. Wang, K. Yasuda, Y. Zhang, S. Liu, K. Watanabe, T. Taniguchi, J. Hone, L. Fu, and P. Jarillo-Herrero, Interfacial ferroelectricity in rhombohedral-stacked bilayer transition metal dichalcogenides, Nature Nanotechnology 17, 367 (2022).
  * Sun  _et al._ [2022] Y. Sun, S. Xu, Z. Xu, J. Tian, M. Bai, Z. Qi, Y. Niu, H. H. Aung, X. Xiong, J. Han, C. Lu, J. Yin, S. Wang, Q. Chen, R. Tenne, A. Zak, and Y. Guo, Mesoscopic sliding ferroelectricity enabled photovoltaic random access memory for material-level artificial vision system, Nature Communications 13, 5391 (2022).
  * Wan  _et al._ [2022] Y. Wan, T. Hu, X. Mao, J. Fu, K. Yuan, Y. Song, X. Gan, X. Xu, M. Xue, X. Cheng, C. Huang, J. Yang, L. Dai, H. Zeng, and E. Kan, Room-temperature ferroelectricity in $1{\mathrm{T}}^{{}^{\prime}}$-${\mathrm{ReS}}_{2}$ multilayers, Phys. Rev. Lett. 128, 067601 (2022).
  * Miao  _et al._ [2022] L.-P. Miao, N. Ding, N. Wang, C. Shi, H.-Y. Ye, L. Li, Y.-F. Yao, S. Dong, and Y. Zhang, Direct observation of geometric and sliding ferroelectricity in an amphidynamic crystal, Nature Materials 21, 1158 (2022).
  * Yasuda  _et al._ [2024] K. Yasuda, E. Zalys-Geller, X. Wang, D. Bennett, S. S. Cheema, K. Watanabe, T. Taniguchi, E. Kaxiras, P. Jarillo-Herrero, and R. Ashoori, Ultrafast high-endurance memory based on sliding ferroelectrics, Science 0, eadp3575 (2024).
  * Bian  _et al._ [2024] R. Bian, R. He, E. Pan, Z. Li, G. Cao, P. Meng, J. Chen, Q. Liu, Z. Zhong, W. Li, and F. Liu, Developing fatigue-resistant ferroelectrics using interlayer sliding switching, Science 0, eado1744 (2024), https://www.science.org/doi/pdf/10.1126/science.ado1744 .
  * Yang  _et al._ [2018] Q. Yang, M. Wu, and J. Li, Origin of two-dimensional vertical ferroelectricity in WTe2 bilayer and multilayer, The Journal of Physical Chemistry Letters 9, 7160 (2018).
  * Wu [2021] M. Wu, Two-dimensional van der Waals ferroelectrics: Scientific and technological opportunities, ACS Nano 15, 9229 (2021).
  * Bennett and Ghosez [2024] D. Bennett and P. Ghosez, Asymmetric dynamical charges in two-dimensional ferroelectrics (2024), arXiv:2404.10549 [cond-mat.mtrl-sci] .
  * He  _et al._ [2024] R. He, B. Zhang, H. Wang, L. Li, P. Tang, G. Bauer, and Z. Zhong, Ultrafast switching dynamics of the ferroelectric order in stacking-engineered ferroelectrics, Acta Materialia 262, 119416 (2024).
  * [34] See supplemental material at ** for Born effective charges of BN bilayer, computational details, model accuracy, analytical calculation of coercive electric field, and supplemental figures, which include refs.[11,24,35-50].
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Physical Review B - Condensed Matter and Materials Physics 54, 11169 (1996), cited by: 84293.
  * Blöchl [1994] P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B 50, 17953 (1994).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Grimme [2006] S. Grimme, Semiempirical GGA-type density functional constructed with a long-range dispersion correction, Journal of Computational Chemistry 27, 1787 (2006), https://onlinelibrary.wiley.com/doi/pdf/10.1002/jcc.20495 .
  * Wang  _et al._ [2019] V. Wang, N. Xu, J. C. Liu, G. Tang, and W. Geng, VASPKIT: A user-friendly interface facilitating high-throughput computing and analysis using vasp code, Comput. Phys. Commun. 267, 108033 (2019).
  * Jinnouchi  _et al._ [2019a] R. Jinnouchi, J. Lahnsteiner, F. Karsai, G. Kresse, and M. Bokdam, Phase transitions of hybrid perovskites simulated by machine-learning force fields trained on the fly with Bayesian inference, Phys. Rev. Lett. 122, 225701 (2019a).
  * Jinnouchi  _et al._ [2019b] R. Jinnouchi, F. Karsai, and G. Kresse, On-the-fly machine learning force field generation: Application to melting points, Phys. Rev. B 100, 014105 (2019b).
  * Jinnouchi  _et al._ [2020] R. Jinnouchi, F. Karsai, C. Verdi, R. Asahi, and G. Kresse, Descriptors representing two- and three-body atomic distributions and their effects on the accuracy of machine-learned inter-atomic potentials, The Journal of Chemical Physics 152, 234102 (2020), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/5.0009491/15575269/234102_1_online.pdf .
  * Allen and Tildesley [2017] M. P. Allen and D. J. Tildesley,  _Computer Simulation of Liquids_ (Oxford University Press, 2017).
  * Gonze and Lee [1997] X. Gonze and C. Lee, Dynamical matrices, Born effective charges, dielectric permittivity tensors, and interatomic force constants from density-functional perturbation theory, Phys. Rev. B 55, 10355 (1997).
  * Musaelian  _et al._ [2023] A. Musaelian, S. Batzner, A. Johansson, L. Sun, C. J. Owen, M. Kornbluth, and B. Kozinsky, Learning local equivariant representations for large-scale atomistic dynamics, Nature Communications 14, 579 (2023).
  * Yu  _et al._ [2024] H. Yu, S. Deng, M. Xie, Y. Zhang, X. Shi, J. Zhong, C. He, and H. Xiang, Switchable ferroelectricity in subnano silicon thin films (2024), arXiv:2407.01914 [cond-mat.mtrl-sci] .
  * Thompson  _et al._ [2022] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, Lammps - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Computer Physics Communications 271, 108171 (2022).
  * Nosé [1984] S. Nosé, A unified formulation of the constant temperature molecular dynamics methods, The Journal of Chemical Physics 81, 511 (1984), cited by: 13724; All Open Access, Green Open Access.
  * OriginLab Corporation [2024] OriginLab Corporation, Origin (2024).
  * Weston  _et al._ [2020] A. Weston, Y. Zou, V. Enaldiev, A. Summerfield, N. Clark, V. Zólyomi, A. Graham, C. Yelgel, S. Magorrian, M. Zhou, J. Zultak, D. Hopkinson, A. Barinov, T. H. Bointon, A. Kretinin, N. R. Wilson, P. H. Beton, V. I. Fal’ko, S. J. Haigh, and R. Gorbachev, Atomic reconstruction in twisted bilayers of transition metal dichalcogenides, Nature Nanotechnology 15, 592 (2020).
  * Shi  _et al._ [2023] B. Shi, H. Wang, W. Jiang, Y. Feng, P. Guo, H. Gao, Z. Gao, and W. Ren, Electric field effect of sliding graphene/hexagonal boron nitride heterobilayer, Applied Surface Science 636, 157816 (2023).
  * King-Smith and Vanderbilt [1993] R. D. King-Smith and D. Vanderbilt, Theory of polarization of crystalline solids, Phys. Rev. B 47, 1651 (1993).
  * Resta [1992] R. Resta, Theory of the electric polarization in crystals, Ferroelectrics 136, 51 (1992), https://doi.org/10.1080/00150199208016065 .
  * Batzner  _et al._ [2022] S. Batzner, A. Musaelian, L. Sun, M. Geiger, J. P. Mailoa, M. Kornbluth, N. Molinari, T. E. Smidt, and B. Kozinsky, E(3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials, Nature Communications 13, 2453 (2022).
  * Gasteiger  _et al._ [2022] J. Gasteiger, J. Groß, and S. Günnemann, Directional message passing for molecular graphs (2022), arXiv:2003.03123 [cs.LG] .
  * Landauer [1957] R. Landauer, Electrostatic considerations in BaTiO3 domain formation during polarization reversal, Journal of Applied Physics 28, 227 (1957), https://pubs.aip.org/aip/jap/article-pdf/28/2/227/18316454/227_1_online.pdf .
  * Xu  _et al._ [2017a] B. Xu, V. Garcia, S. Fusil, M. Bibes, and L. Bellaiche, Intrinsic polarization switching mechanisms in ${\mathrm{BiFeO}}_{3}$, Phys. Rev. B 95, 104104 (2017a).
  * Xu  _et al._ [2017b] B. Xu, J. Íñiguez, and L. Bellaiche, Designing lead-free antiferroelectrics for energy storage, Nature Communications 8, 15682 (2017b).
  * Daumont  _et al._ [2012] C. Daumont, W. Ren, I. C. Infante, S. Lisenkov, J. Allibe, C. Carrétéro, S. Fusil, E. Jacquet, T. Bouvet, F. Bouamrane, S. Prosandeev, G. Geneste, B. Dkhil, L. Bellaiche, A. Barthélémy, and M. Bibes, Strain dependence of polarization and piezoelectric response in epitaxial BiFeO3 thin films, Journal of Physics: Condensed Matter 24, 162202 (2012).
  * Lisenkov  _et al._ [2009] S. Lisenkov, D. Rahmedov, and L. Bellaiche, Electric-field-induced paths in multiferroic ${\mathrm{BiFeO}}_{3}$ from atomistic simulations, Phys. Rev. Lett. 103, 047204 (2009).

## SI. Born effective charges (BECs) of BN bilayer

The force acting on ion $i$ in the presence of a finite electric field $E$ can
be determined using the formula $\bm{F}_{i}=Z_{i}(\bm{R})\bm{E}$, where
$Z_{i}(\bm{R})$ denotes the zero-field Born effective change (BEC) tensor of
ion $i$ at the coordinates $\bm{R}$. During polarization switching in sliding
ferroelectrics, that the positional rearrangement among ions within the same
layer remains negligible, with the primary change attributed to interlayer
sliding. Consequently, we group two adjacent $\mathrm{B}$ and $\mathrm{N}$
ions within the same layer as a collective entity to examine the BECs
distribution along the sliding pathway. Illustratively, as depicted in Fig.
S1, ions experience forces of equal magnitude but opposite direction under the
influence of an external field. Specifically, when $E_{\perp}$ is applied, the
distribution of forces is symmetric with respect to the transition state SP
[see Fig. S1(a)]. Conversely, with $E_{parallel}$ applied in the sliding
direction, the distribution becomes antisymmetric [see Fig. S1(b)]. In both
situations, the BECs exhibit significant variations along the sliding path
from BA to AB. Previous MD simulations that use fixed BEC values cannot
precisely capture the behavior of the $\mathrm{BN}$ bilayer under these
electric field conditions. This limitation is particularly evident for
$E_{parallel}$, where it cannot calculate the force by employing the average
BEC value along the sliding pathway, which would be zero. Hence, our proposed
model, capable of predicting BECs for varying coordinates $\bm{R}$, proves to
be indispensable.

## SII. Computational details and model accuracy

We employ the Vienna _ab initio_ Simulation Package (VASP) [35] to perform DFT
calculations using the projector augmented wave (PAW) method [36] and the
generalized gradient approximation (GGA) with the Perdew-Burke-Ernzerh (PBE)
exchange-correlation functional [37]. The plane-wave cutoff energy is set to
520 eV. The cell constant along the $c$-axis is fixed at 30 Å, ensuring a
sufficient vacuum space. A $12\times 12\times 1$ $k$-point mesh and vdW
corrections with the DFT-D2 method of Grimme [38] are applied for structure
optimization. The optimized lattice constant of the AB-stacked $\mathrm{BN}$
bilayer is $a=b=2.511$ Å with an interlayer distance of $3.074$ Å, which is
consistent with previous theoretical calculations [11]. For other parallel-
stacked configurations corresponding to different sliding vectors, the lattice
constant remains almost unchanged, and the interlayer distance increases
slightly (not exceeding 0.35 Å). The primary cells are expanded to $5\times
5\times 1$ supercells using Vaspkit [39], and MD calculations are performed
using on-the-fly machine learning force fields (MLFF) [40, 41, 42] for more
efficient sampling. We use a $3\times 3\times 1$ $k$-point mesh and the NVT
ensemble with a Langevin thermostat [43]. The BECs of each atom in different
configurations are further calculated using density functional perturbation
theory (DFPT) [44].

In our training, the cutoff radius is set to $R_{c}=7\text{\r{A}}$,
encompassing nearly 100 pairs of each atom formed within the 19th nearest
neighbor. The learning rate is designed to decay exponentially from $5\times
10^{-3}$ to $9.7\times 10^{-6}$. For the Allegro [45] model, trained on a
dataset containing 1659 structures from MD sampling blow 3000K, we select a
(64, 128, 256, 256) dimensional embedding network, a (256, 256, 256)
dimensional fitting network, and a (128) dimensional per-edge final network.
This configuration yields mean absolute errors (MAEs) of 0.053 meV/atom for
energy, 5.311 meV/Å for forces and 0.063 meV/Å3 for stress. For our DREAM-
Allegro [46] model (DA), trained on a dataset containing 3484 structures from
MD sampling blow 500K, the loss weight coefficients for energy, forces,
stress, and BEC0 (which measures the system’s electrical neutrality) are all
set to 1, and the coefficient for BEC is increased to 10 due to the relatively
small absolute values of BECs. The per-edge final network dimensions are
adjusted to (128, 128, 64) dimensions, resulting in a model with MAEs of 0.030
meV/atom for energy, 3.851 meV/Å for forces, 0.076 meV/Å3 for stress and 0.008
$|$e$|$ for BEC (after correcting with BEC0). This indicates that our models
are highly accurate. Fig. S2 presents the parity plot of energy and BEC.

The MD simulations are performed using the LAMMPS software package [47] with a
newly integrated atom type named “BEC”. Periodic boundary conditions are
applied along the $a$ and $b$ axes, while a non-periodic, fixed boundary
condition is used along the $c$ axis to accurately represent 2D systems. The
canonical NVT ensemble, with temperature regulation provided by a Nose-Hoover
thermostat [48], is employed. The simulation box is a $15\times 15\times 1$
supercell, consisting of 900 atoms, with an integration time step of 1
femtosecond (fs).

## SIII. Analytical calculation of coercive electric field

The fitting of energies and dipoles per unit cell along the BA to AB sliding
pathway is performed using Origin [49]. Considering their even and odd
function characteristics, $\varepsilon_{0}(t)$ and $p_{\parallel}(t)$ are
fitted with quartic polynomials, and $p_{\perp}(t)$ is fitted with cubic
polynomials.Constraints are applied to ensure that $\varepsilon_{0}$ is
minimized at the BA (AB) state, $p_{\parallel}$ is zero at the BA (AB) state,
and $p_{\perp}$ reaches a maximum (minimum) at the BA (AB) state. The fitting
results are presented in Fig.S3.

When the critical electric field is defined as the minimum field at which the
local energy minimum of $\varepsilon(t)$ disappears, i.e.,
$\mathrm{d}\varepsilon(t)/\mathrm{d}t\leq 0$ along the BA to SP sliding
pathway, the exact analytical expression for $E_{\perp,c}$ is given by:

$\displaystyle
E_{\perp,c}=\frac{-2e_{4}t_{\mathrm{BA}}^{2}-p_{2}^{\parallel}E_{\parallel}}{3t_{\mathrm{BA}}p_{3}^{\perp}}\sqrt{\frac{8x}{m(m-2)}},$
(S6) $\displaystyle\text{where}\
m=\frac{3}{4}+\frac{3}{2}x+\frac{1}{4}\sqrt{36x^{2}-60x+9}$
$\displaystyle\text{and}\
x=\frac{3(e_{4}-p_{4}^{\parallel}E_{\parallel})t_{\mathrm{BA}}^{2}}{2e_{4}t_{\mathrm{BA}}^{2}+p_{2}^{\parallel}E_{\parallel}}.$

Eq.4 approximates Eq. S6 to the order of 0.5 of $E_{\parallel}$.

## Supplemental figures

Figure S1:  The sum of sliding-related BEC components of adjacent $\mathrm{B}$
and $\mathrm{N}$ ions in the top layer (orange line) and bottom layer (blue
line) along the sliding pathway, calculated using DFPT. (a) The BEC component
related to $E_{\perp}$ and the in-plane force in the sliding direction
$F_{\parallel}$ ($F_{\parallel}=Z_{\rm\parallel,\perp}E_{\perp}$). (b) The BEC
component related to $E_{\parallel}$ and $F_{\parallel}$
($F_{\parallel}=Z_{\parallel,\parallel}E_{\parallel}$).

Figure S2:  The parity plot of the DREAM-Allegro (DA) model compared to DFT
results. Comparison of energies $\varepsilon_{0}$ (a) and BECs $Z$ (b) from
the DA model against DFT calculations for all configurations in the training
dataset. The insets display the distribution of the MAEs.

Figure S3:  Fit plots and goodness of fit ($R^{2}$) for energies
$\varepsilon_{0}(t)$ (a), in-plane dipoles $p_{\parallel}(t)$ (b) and out-of-
plane dipoles $p_{\perp}(t)$ (c) per unit cell.

Figure S4:  The variation of the MD simulation obtained $E_{\perp,c}$ with
$E_{\parallel}$ orienting from $0^{\circ}$ to $60^{\circ}$ at a fixed
magnitude of 0.1V/Å at 0.1K. $E_{\perp,c}$ is determined by identifying the SP
state during the sliding process when $E_{\perp}$ along the $-z$ direction is
incrementally increasing over time at a rate of 0.02 V/(Å$\cdot$ps). Within
the range of $\left(0^{\circ},60^{\circ}\right]$, the sliding consistently
occurs at $60^{\circ}$, while at $0^{\circ}$, the sliding occurs randomly at
$60^{\circ}$ or $300^{\circ}$. The closer the direction of $E_{\parallel}$ is
to the direction of sliding, the smaller the $E_{\perp,c}$.

Figure S5:  The variation of the MD simulation obtained $E_{\perp,c}$ with
$E_{\parallel}$ at different magnitudes oriented at $60^{\circ}$ and at
different temperatures. $E_{\perp,c}$ is determined as the minimum vertical
field, adjusted in increments of 0.05 V/Å, that can initiate sliding within 30
ps. The pentagram marks indicate the the minimum $E_{t,c}$ at each
temperature, corresponding to the blue line in Fig.4(c).

Figure S6:  (a) Energy $\varepsilon_{0}$ and polarization $P$ distribution of
$\mathrm{MoS_{2}}$, calculated by the DFT and Berry phase approach. Positions
within the hexagon correspond to different sliding vectors, with the center
representing the ferroelectric XM state [24, 50]. The shades of grey in the
background indicate energy levels. The direction of colored arrows shows in-
plane polarizations $P_{\parallel}$, while their colors represent out-of-plane
polarization $P_{\perp}$. Energy $\varepsilon_{0}$ (b) and polarization $P$
(c) along the designed sliding pathway.