Deterministic and Efficient Switching of Sliding Ferroelectrics

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.

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].

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

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

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.

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}$.

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/Å⁴, $e_{0}=7.41$ meV, $p_{4}^{\parallel}=5.82\times 10^{-3}$ e/ų, $p_{2}^{\parallel}=-1.66\times 10^{-2}$ e/Å, and $p_{3}^{\perp}=9.13\times 10^{-3}$ e/Ų [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

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

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

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.