Ferroelectricity in van der Waals multilayers via interfacial polarization engineering

Ferroelectrics, characterized by two switchable polarization states representing logic “0” and “1”, stand out as ideal materials for developing non-volatile random-access memories^1,2. Over the past decade, the growing demand for device miniaturization has aroused tremendous interest in studying low-dimensional van der Waals (vdW) ferroelectric (FE) materials, including In₂Se₃ and analogous III₂-VI₃ materials^3,4, CuInP₂S₆⁵, 1T-MoTe₂⁶, GaSe⁷ and group-IV monochalcogenide monolayers⁸. However, compared to the vast number of reported two-dimensional (2D) materials, those with intrinsic ferroelectricity remain scarce. In 2017, Li and Wu proposed the concept of sliding ferroelectricity⁹, which provides a new paradigm to construct 2D ferroelectrics using non-polarized parent materials. The polarization in sliding ferroelectrics is built from the net interlayer electric dipoles and can be switched via the relative sliding of stacked layers^10,11. Later, the sliding ferroelectricity has been experimentally characterized in various materials, including parallel or antiparallel BN^12,13,14, transition-metal dichalcogenides (TMDs)^15,16, γ-InSe¹⁷, as well as multiwalled nanotubes¹⁸ and 0D nanotube crossbars¹⁹. The sliding ferroelectrics are also distinguished as a unique platform for exploring the interplay between ferroelectricity and other physical properties, such as ferromagnetism²⁰, ferrovalley²¹, moiré excitons²², nonlinear optical frequency conversion²³, bulk photovoltaic effects^24,25 and layer Hall effect^26,27, which greatly enriches the fundamental interest and potential applications on this topic^28,29,30,31.

As the sliding ferroelectricity is extensively studied, ferroelectric vdW multilayers are recognized as a host to even more interesting behaviors. The intricate interplay among stacking number, order, and atomic species in vdW multilayers can introduce peculiar ferroelectricity, such as tunable multiple polarization states^32,33,34,35 and atypical ferroelectricity induced by the symmetry breaking across layers^36,37. Recently, intrinsically nonpolar Bernal-stacked bilayer graphene (BLG) exhibits ferroelectricity when sandwiched between hexagonal boron nitride (hBN), a phenomenon attributed to moiré potential of the vdW heterostructures and electron correlation effects in the graphene (Gr) bilayer^38,39,40. The ferroelectric hysteresis of these heterostructures can be sustained even at room temperature. A subsequent theory suggests a concept of across-layer asymmetric sliding ferroelectricity³⁶, but the estimated polarization still does not align with the experimental results. Given the emergence of compelling phenomena in the hBN-capsulated BLG⁴¹, such as non-volatile ratcheting electronic states, it is imperative to thoroughly elucidate the origin of robust electric polarizations in vdW multilayers. Furthermore, understanding how each layer influences interfacial charge transfer and polarization in adjacent layers is essential for constructing FE multilayers with customizable polarization; yet this issue remains unaddressed.

Here, we report an intrinsic effect of polarization accumulation in vdW multilayers utilizing Gr-capsulated BN layers as an example. We first demonstrate this effect by showing an unusual FE phase in a seemingly non-polar system comprised of a BN monolayer sandwiched between Gr monolayers. The overall polarization increases with inserting more layers of rhombohedral BN (rBN) between the Gr layers, exhibiting a cumulative effect. Charge density analysis shows that intralayer electric dipoles of component layers are enhanced by interlayer charge redistribution, accompanied by cross-layer electrostatic interaction. Such an interaction reduces the interlayer spacing and enhances the interlayer polarization. Thickening the Gr layers can further enhance interfacial polarization due to decreased spacing at the Gr/BN interfaces and the additionally induced polarization within the Gr layers. The pathway for ferroelectric switching is analyzed based on minimizing the generalized stacking fault energy (GSFE) barrier, which manifests as sequential or cooperative interlayer slipping driven by an external field. The ferroelectric polarization within the Gr-capsulated BN layers exerts a significant modulation effect on the doping level of the surface Gr layers, thereby enabling the realization of multilevel conductive states. These results can be extended to other multilayer systems, such as incommensurate molybdenum disulfide MoS₂/rBN/MoS₂ structures, opening a new avenue for designing multilayer FE devices through cross-layer interactions.

Ferroelectricity in graphene-encapsulated boron nitride

We first focus on a trilayer system comprising a BN monolayer sandwiched between two Gr monolayers, denoted as Gr/BN/Gr. In the ground state of the Gr/BN bilayer³⁶, boron (B) atoms are located above carbon (C) atoms while nitrogen (N) atoms sit atop the centers of carbon hexagons. We refer to this interface as C-B for convenience. The ground state of Gr/BN/Gr features two identical C-B interfaces, arranged in an ABA stacking order with mirror symmetry, resulting in no out-of-plane polarization, as illustrated in Supplementary Fig. 1. Hereafter, polarization refers specifically to out-of-plane polarization unless otherwise noted. Different stacking phases of AB, BA, and AA coexist at an incommensurate Gr/BN interface⁴², and by slipping one of the Gr layers in the ABA configuration, the system transitions into an electrically polarized ABC mode due to symmetry breaking (Fig. 1a). At the new Gr/BN interface, the N atoms are atop the C atoms while the B atoms are situated above the centers of carbon hexagons, denoted as the N-C interface. The FE polarization of the ABC phase can be switched to the CBA phase through the cooperative slipping of both the top and bottom Gr layers (indicated by red arrows in Fig. 1a).

a Schematic of Gr/BN/Gr with ABC and CBA stacking phases. FE polarizations (gray arrows) are induced within the asymmetric phases, as reflected by the vertical deviation of the BN layer from the central plane of the structure. The polarization can be reversed by cooperative slipping of the top and bottom Gr layers (red arrows). Interlayer spacings (d₁ and d₂) and corresponding van der Waals (vdW) energies (E₁ and E₂) are denoted for the upper and lower Gr/BN interfaces. b Differences in d₁ and d₂ between ABA and ABC phases. c vdW energy and interlayer spacing as functions of normalized transformation coordinate.

Full size image

Interfacial sliding not only changes the stacking order but also alters the interlayer spacing in the Gr/BN/Gr system. Here, we denote the interlayer spacing between BN and the top Gr layer as d₁, and that between BN and the bottom Gr layer as d₂. As the stacking order shifts from ABA to ABC, d₂ increases from 3.32 to 3.48 Å, while d₁ remains nearly unchanged at 3.32 Å (Fig. 1b). In the ABC configuration, d₁ (i.e., C-B spacing) and d₂ (i.e., N-C spacing) agree with previously calculated results of 3.3 and 3.5 Å using the accurate adiabatic fluctuation-dissipation theorem method⁴³. These values also align with the experimentally measured 3.32 ± 0.43 Å for the hBN/Gr interface⁴⁴. Different d₁ and d₂ indicate that the sandwiched BN layer deviates vertically from the central plane of the polarized ABC phase. This deviation can be understood from the interfacial vdW energy, defined as ({E}_{{rm{vdW}}}={E}_{{rm{total}}}-{E}_{{rm{a}}}-{E}_{{rm{b}}}), where ({E}_{{rm{vdW}}}) is the total energy of the multilayer, and ({E}_{{rm{a}}}) and ({E}_{{rm{b}}}) are the energies of the individual layers above and below the interface, respectively. The calculated vdW energy profiles for both the upper and lower BN/Gr interfaces along the polarization switching pathway closely reproduce the variation trends of the interlayer spacings d₁ and d₂ (Fig. 1c).

The key parameters characterizing sliding ferroelectricity are the ferroelectric polarization P and the switching barrier. P of the Gr-sandwiched BN monolayer calculated using the Berry-phase method⁴⁵ with dipole correction⁴⁶ is 1.01 pC/m (Fig. 2a), somewhat surprising since a BN monolayer lacks an out-of-plane polarity. The calculated P is even higher than the value of 0.91 pC/m for an rBN bilayer. The energy barrier at the saddle point (SP) is 4.4 meV per unit cell (u.c), slightly lower than 6.4 meV/u.c for the across-layer ferroelectricity in the previously studied BN/Gr/Gr/BN system³⁶. Interestingly, defining the vertical deviation of the BN layer as (triangle d=({d}_{1}-{d}_{2})/)2(,) we obtain a linear relationship between Δd and P along the FE switching pathway (Supplementary Fig. 2).

a Energy and polarization as functions of normalized transformation coordinate, where the structures of initial (ABC), saddle-point (SP), and final (CBA) phases are illustrated in the insets. b Isosurface plots (1.2 × 10⁻⁵ e/ų) of differential charge density distribution of the ABC stacking. Electron accumulation and depletion regions are colored yellow and blue, respectively.

Full size image

It is known that charge transfer occurs in a heterostructure to balance the Fermi level at the interface. To unravel the origin of polarization in the ABC-stacked Gr/BN/Gr, we calculate the charge density difference by subtracting the sum of the individual monolayer charge densities from that of the entire system (Fig. 2b). At the upper Gr/BN interface, electrons accumulate between the vertically aligned B and C atoms but deplete near the carbon hexagons. This distinct vertical separation of charges creates a downward-pointing dipole moment (or ferroelectric field), corresponding to an upward polarization. At the lower Gr/BN interface, electrons are depleted near the C atoms atop N atoms but accumulate near the C atoms of the other sublattice. Yet, the vertical charge transfer is less pronounced at this interface than at the upper one, indicating a weaker interfacial polarization. Our Bader charge analysis (Supplementary Fig. 3) verifies that during the transition from ABC to SP stacking, the top and bottom Gr layers lose and gain electrons, respectively. Intriguingly, unlike the screening effect normally observed in Gr multilayers⁴⁷, the Gr monolayers at both interfaces are critical for building up the ferroelectric field. Band structures indicate that the two Gr layers are oppositely doped in the ABC stacking (Supplementary Fig. 4). As the overall polarization gradually diminishes following the transition to the SP phase, the band structures of the two Gr layers converge to the same.

Cumulative polarization and electrostatic coupling

The notable ferroelectricity in the ABC-stacked Gr/BN/Gr structure prompts us to investigate how interfacial polarization is influenced by BN thickness. The sandwiched BN monolayer in Fig. 1a is replaced with n layers of rBN (n-rBN), while preserving the C-B and N-C interfaces, thus forming a Gr/n-rBN/Gr structure (Fig. 3a). The calculated polarization (P) values of both the Gr/n-rBN/Gr and bare n-rBN systems are plotted in Fig. 3b. As n increases from 1 to 4, the P values for bare n-rBN are 0, 0.91, 1.94, and 3.00 pC/m, while those for Gr/n-rBN/Gr rise to 1.01, 2.05, 3.10, and 4.16 pC/m, respectively. Both bare n-rBN and Gr/n-rBN/Gr exhibit a nearly linear increase in polarization with the addition of BN layers, highlighting the cumulative effect of interfacial ferroelectricity. Notably, the polarization of the Gr/2-rBN/Gr shows a 124% enhancement compared to that of the bare 2-rBN, a level of improvement that holds true across various vdW correction methods (Supplementary Fig. 5). Moreover, the polarization difference (ΔP) between n-rBN and Gr/n-rBN/Gr rises as n increases from 1 to 3. This indicates that ΔP arises not only from the newly formed polarizations at the Gr/BN interfaces, but also from additional effects within the sandwiched rBN.

a Schematic of the Gr/n-rBN/Gr multilayer stacking, with the gray arrow indicating the polarization. b Polarizations of Gr/n-rBN/Gr and bare n-rBN (without Gr), along with their difference ΔP. c Comparison of the N-B spacings, where d_xy denotes the interlayer spacing between BN layers x and y. d Negative piezoelectric responses under applied strains for bilayers with N-B, C-B, and N-C interfaces (insets). Dashed lines on the strain-polarization curves indicate the maximum strains for different interfaces in Gr/n-rBN/Gr (n ≤ 4). e Overall polarization of Gr/n-rBN/Gr with the cumulative interfacial polarizations obtained through three approaches. f Band structures of Gr/n-rBN/Gr near the Fermi level.

Full size image

To clarify the variation in ΔP, we investigated the interlayer spacings within n-rBN and Gr/n-rBN/Gr. Figure 3c shows the interlayer spacings between the BN layers (N-B spacings) in both systems, denoted as d_xy, where x and y correspond to the indices of the BN monolayers shown in Fig. 3a. In 2-rBN, the value of d₁₂ is 3.342 Å, while in 3-rBN, both d₁₂ and d₂₃ shrink by approximately 0.015 Å, reaching 3.328 and 3.326 Å, respectively. The spacing compressions suggest an attractive force between the 1^st and 3^rd BN layers in 3-rBN, induced by the second neighboring BN layer spanning the central layer. As n increases to 4, d₁₂, d₂₃, and d₃₄ are 3.325, 3.310, and 3.323 Å, respectively. Notably, d₂₃ is shorter than d₁₂ and d₃₄ in the 4-rBN structure, as the central interface becomes compressed by the second neighboring layers on both sides. The outer interlayer spacings in 4-rBN (d₁₂ and d₃₄) are only 0.003 Å smaller than the corresponding d₁₂ and d₂₃ in 3-rBN, suggesting a weaker influence from the third neighboring layer. With the n-rBN coated by two heterogeneous Gr monolayers, all the N-B spacings are uniformly compressed to about 3.30 Å. Furthermore, the interlayer spacing at the central interface of the Gr/2-rBN/Gr (3.302 Å) system is shorter than that of the 4-rBN system (3.310 Å), suggesting that the Gr monolayers facilitate stronger electrostatic interactions across layers. We also examine the interlayer spacings of C-B and N-C interfaces in the Gr/n-rBN/Gr multilayer and Gr/BN bilayer (Supplementary Fig. 6). The spacings in the multilayers are reduced compared to the bilayer, indicating a compressive effect at the heterogeneous interfaces.

Compression of a polarized interface can enhance its polarization, known as the negative piezoelectric effect⁴⁸. To quantitatively evaluate this effect in the Gr/n-rBN/Gr multilayer, we divide the system into three components: the N-B, C-B, and N-C interfaces. The strain-induced polarization variations for these interfaces are illustrated in Fig. 3d. At zero strain, the polarizations for the N-B, C-B, and N-C interfaces are 0.91, 0.89, and 0.08 pC/m, respectively. The corresponding values increase linearly to 1.12, 1.04, and 0.16 pC/m as the strain is reduced to -1.5%. We then compare the overall polarization of the Gr/n-rBN/Gr system with the summed polarization of its component interfaces, as shown in Fig. 3e. The polarizations for the component interfaces are computed using three methods: (1) the polarization of the unstressed interface, (2) the polarization considering interfacial strain and the corresponding strain-polarization relationship from Fig. 3d, and (3) the polarization using fixed atomic positions from the multilayer. The results show that the overall polarization aligns well with the accumulation of polarizations at all interfaces, incorporating the piezoelectric effect, though the interfacial strains are induced by cross-layer interactions rather than by external stress. Specifically, the N-B spacing in Gr/2-rBN/Gr is 3.300 Å, corresponding to a strain of −1.26%. This strain induces a polarization enhancement of 0.17 pC/m, accounting for an 18% increase relative to the interfacial polarization of bare 2-rBN (0.91 pC/m). The piezoelectric effect, in tandem with the newly formed polarization at the Gr/BN interfaces, pinpoints the ΔP shown in Fig. 3b. We further investigate the impact of the number of BN layers on Gr doping in the Gr/n-rBN/Gr structure, with the band structures for increasing n shown in Fig. 3f. It is observed that the Dirac points of both Gr layers transform into Dirac loops. The radii of these Dirac loops increase as n increases, reaching a maximum value around n = 3. Since the interlayer spacings and polarizations at both Gr/BN interfaces remain constant with increasing n, the enhanced doping in Gr highlights the dominant role of cross-layer electrostatic interactions.

As the number of BN layers in n-rBN and Gr/n-rBN/Gr increases further (n ≥ 5), the minimum N-B spacings in both systems progressively converge towards the bulk rBN interlayer spacing of 3.27 Å (Supplementary Fig. 7a). The polarization of n-rBN shows a linear increase with n (Supplementary Fig. 7b), consistent with the cumulative law of interfacial polarization. Similarly, the overall polarization of the Gr/n-rBN/Gr system is projected to continuously increase with n, though the system transitions to a metallic state for n ≥ 5, which prevents polarization calculations using the Berry-phase method. It is important to note, however, that the surface properties of a ferroelectric multilayer are expected to saturate at a certain thickness. Recently, Cao et al. experimentally demonstrated that the surface potential of rhombohedral MoS₂ increases linearly with thickness up to 7 layers and saturates around 10 layers⁴⁹. They further predicted that the surface potential of rBN would saturate at approximately 40 layers, attributed to polarization-induced bandgap closure that leads to the appearance of free carriers. Based on these thickness-dependent variations in surface potential and our observations of Gr doping, we suggest that until the system undergoes a metallic phase transition, the surface properties of a multilayer ferroelectric exhibit a direct dependence on the overall polarization.

Now we elucidate the mechanism underlying the compressed spacing and enhanced polarization at the interfaces within vdW multilayers, with the rBN system serving as an illustrative example. Despite the overall non-polarity of a BN monolayer, the electronegativity difference between B and N atoms leads to the formation of local electric dipoles (N-B dipoles), as illustrated in Fig. 4a. When two BN layers are stacked together to form an N-B interface, electron redistribution occurs at the interface, which enhances the intralayer N-B dipoles in both layers compared to a single BN layer (Fig. 4b). Hirshfeld-I charge^50,51 analysis (Supplementary Fig. 8) supports this viewpoint. In an isolated BN monolayer, B and N atoms carry charges of 1.266 e and −1.266 e, respectively. However, in an rBN bilayer, the absolute values of the charges on B and N atoms increase to about 1.291 e. When a third BN layer is stacked onto the bilayer to form a 3-rBN system, the enhanced intralayer N-B dipoles in the original bilayer promote charge transfer with the newly stacked BN layer, demonstrating a cross-layer electrostatic interaction (Fig. 4c). The enhanced interlayer charge transfer in the 3-rBN system reduces spacings and enhances polarizations at both N-B interfaces. As the number of layers increases to four (Fig. 4d), the central interface, encapsulated by second neighboring layers on both sides, undergoes a doubled amplification of interlayer compression and polarization enhancement. This is evidenced by the B and N atomic charges at the central interface, which reach 1.313 e and −1.314 e, respectively. As such, interlayer charge transfer and intralayer local dipoles mutually enhance each other, leading to a cross-layer electrostatic coupling effect. Similar interlayer compression is obtained in rhombohedral MoS₂ (Supplementary Fig. 9) and even incommensurate MoS₂/2-rBN/MoS₂ systems (Supplementary Fig. 10), indicating that cross-layer electrostatic coupling is a universal phenomenon in ferroelectric multilayers.

a Local intralayer electric dipoles within a BN monolayer. Electrostatic interactions and interlayer polarizations in (b) bilayer, (c) trilayer, and (d) quadrilayer rBN structures. Red and blue regions denote negative and positive charges of the local N-B dipoles, respectively, while red and orange arrows represent electron transfer and FE polarization. NLs refer to neighboring layers, and 2^nd NL indicates the second neighboring layer.

Full size image

Interfacial polarization involving multilayer graphene

Inspired by the significant ferroelectricity experimentally observed in the hBN-capsulated BLG^38,39, we investigate whether thickening the Gr layers can enhance the polarization of Gr/n-rBN/Gr. As illustrated in Fig. 5a, m₁ and m₂ layers of rhombohedral Gr (rGr) are stacked adjacent to the C-B and N-C interfaces, respectively. The heterostructure is arranged in an ABC stacking order and denoted as m₁/n/m₂ for convenience. Considering the distinct charge transfer behaviors at the C-B and N-C interfaces, we calculate the polarizations for the m/n/1 and 1/n/m systems and analyze the differences in their polarizations compared to the 1/1/1 system (ΔP), as shown in Fig. 5b. The results show that as the number of Gr layers (m) increases, m/n/1 exhibits significant polarization enhancements, while 1/n/m shows negligible changes. Specifically, the ΔP values for m/1/1, m/2/1, and m/3/1 at m = 1 are 1.01, 1.13, and 1.16 pC/m, respectively. These values rise rapidly to 1.19, 1.32, and 1.36 pC/m at m = 2, and further to 1.24, 1.38, and 1.41 pC/m at m = 3. For m > 3, both the m/3/1 and 1/3/m systems transition into metallic states, rendering polarization calculation invalid using the Berry-phase method. We then examine the variations of the C-B spacing in m/n/1 and the N-C spacing in 1/n/m (Fig. 5c), which decrease obviously as m increases, indicating the cross-layer compressive effect from the multilayer Gr. When replacing the sandwiched rBN in m₁/n/m₂ with hBN while maintaining the corresponding C-B and N-C interfaces, similar trends in spacing and polarization are observed (Supplementary Fig. 11), confirming that the ΔP originates from the Gr/BN or Gr/Gr interfaces.

a Schematic of the multilayer stacking, with the gray arrow indicating the polarization. b Polarization differences (ΔP) between m/n/1 and n-rBN (solid squares), and between 1/n/m and n-rBN (hollow squares). c C-B and N-C spacings in the n/m/1 and 1/m/n systems, respectively. d Polarization reduction of the non-3R-like (marked in red) stacks compared to the 3R-like stacks. e Schematic illustration of the dipole model in Gr/2-rBN/Gr, with black arrows indicating the vertically aligned C-C atom pairs.

Full size image

To determine whether polarization exists within the Gr multilayer, we analyze the composition of ΔP (0.18 pC/m) in the 2/1/1 system. The primary piezoelectric effect occurs at the C-B interface, contributing 0.09 pC/m to the ΔP. According to the cumulative law of interfacial polarization, the remaining 0.09 pC/m is thought to originate from the Gr/Gr interface. Additionally, the stacking order of the Gr multilayer affects polarization. Defining a heterogeneous multilayer arranged in ABC order as a 3R-like stack, we find that 3R-like stacks exhibit larger polarizations than their non-3R-like counterparts (Fig. 5d). These differences arise from the characteristics of charge distribution in Gr multilayers, where the p_z orbitals of vertically aligned carbon-carbon (C-C) atom pairs repel the transferred electrons (Supplementary Figs. 12, 13). We can also qualitatively compare the polarizations associated with different Gr stacking orders by applying the dipole model. For instance, in the 2/1/1 system with CABC stacking (Fig. 5e), the transferred electrons in the bilayer Gr are repelled toward the C atoms directly above B atoms. Consequently, the 3R-like stacking order promotes interlayer charge transfer at the Gr/BN interface, while sustaining the interlayer dipoles within the bilayer Gr. In contrast, in BABC (non-3R-like) stacking, the repelled electrons are not aligned above the B atoms, and the vertical alignment of the C-C atom pairs with B atoms diminishes the dielectric performance of the bilayer Gr, leading to a decrease in polarization. Importantly, as the unconventional ferroelectricity measured in hBN-encapsulated BLG has not been well explained, we propose that the interfacial polarization of a 3R-like BN/BLG/BN system can reach 1.31 pC/m or 0.13 μC/cm² (Supplementary Fig. 14), which is comparable to the experimental value of 0.18 μC/cm².

Multilayer polarization switching and Fermi level tuning

Previous experimental studies have demonstrated cumulative polarization and ladder ferroelectricity in TMDs³⁵, and attempted to clarify related atomic-level interlayer sliding in multilayer γ-InSe⁵². However, the polarization switching pathways were suggested from first-principles calculations without a well-defined criterion, hindering precise design and control of interface-engineered ferroelectric devices. Here, we demonstrate a complete ferroelectric reversal mechanism in Gr/2-rBN/Gr and seek to propose general guidelines for interfacial polarization switching in vdW multilayers. For the Gr/2-rBN/Gr multilayer, the full ferroelectric reversal involves slipping all three interfaces: N-C to C-B, N-B reversal, and C-B to N-C. Supplementary Fig. 15 illustrates the generalized stacking fault energy (GSFE) in the multilayer for sequential switching (one interface at a time). It shows that the stacking energy decreases by 15.1 meV/u.c. during the N-C to C-B transition (and vice versa), while the energy barrier for the N-B reversal is 3.6 meV/u.c. Considering the corresponding interlayer slipping in individual Gr/BN and BN/BN bilayers, the energies change moderately to 14.2 and 3.3 meV/u.c., respectively, consistent with previously reported values^43,53. The small GSFE discrepancy between the multilayer (E) and uncoupled bilayers (E_u) suggests that the interlayer electrostatic interactions have a minimal impact on slipping barriers. Therefore, the slipping pathways in multilayers can be effectively decomposed into those of individual bilayers or their combinations.

Compared with the sequential slipping of individual Gr/BN interfaces in the Gr/2-rBN/Gr multilayer, cooperative slipping of both Gr/BN interfaces results in a lower energy barrier of 4.4 meV/u.c., offering a more energetically favorable switching mode. Thus, the full ferroelectric reversal can proceed in two steps (Fig. 6a, b). First, slipping occurs at the BN/BN interface, transitioning the stacking from ABCA to ABAB and reducing the overall polarization from 2.08 pC/m to 0.05 pC/m. Next, both Gr/BN interfaces slip cooperatively, transitioning the stacking from ABAB to CBAC and fully reversing the polarization to −2.08 pC/m. Since the GSFEs of multiple interfaces are nearly decoupled, we can establish guidelines for determining ferroelectric switching pathways in multilayers by minimizing the energy barriers. For homogeneous multilayers, the GSFE is uniform at each interface, resulting in a lower energy barrier for sequential slipping than for cooperative slipping. However, for heterogeneous multilayers, the varying sliding energy profiles across interfaces may favor cooperative slipping to further reduce switching barriers. Notably, applying an external field to an interface alters the GSFE. For instance, in the Gr/BN bilayer, an electric field of 3.5 V/Å enables the C-B interface to overcome the GSFE barrier and transition to the N-C interface (Supplementary Fig. 16). With the cooperative slipping of another Gr/BN (N-C) interface in the Gr/BN/Gr system, the C-B interface needs a lower electric field of 2.5 V/Å for polarization switching. One should note that sequential and cooperative slipping represent ideal switching pathways for ferroelectric multilayers. In practical scenarios, however, the subtle slipping processes associated with polarization changes are susceptible to factors such as the elimination of stacking faults¹⁷, strain-induced out-of-plane distortions⁵⁴, and even vacancy defects⁵⁵. These factors, although increasing the complexity of the switching mechanism, also enrich the system by introducing diverse pinning strategies, thereby enabling the realization of multiple metastable polarization states.

a Structures of initial (ABCA), middle (ABAB) and final (CBAC) stacking orders during switching. The red and orange arrows indicate interlayer sliding and polarizations, respectively. b Energy and polarization as functions of the normalized transformation coordinate, with SP1 and SP2 denoting two saddle points, and ΔE_u representing the uncoupled energy. c Band structures of typical configurations along the polarization switching pathway.

Full size image

During the ferroelectric reversal process of Gr/2-rBN/Gr, the overall polarization decreases monotonically, leading to the emergence of multi-polarized states (Fig. 6b) with effective pinning effects from moiré interfaces^16,56,57. The evolution of band structures along the switching pathway is shown in Fig. 6c. Despite sharing the C-B and N-C interfaces, the ABCA, the first saddle point (SP1), and the ABAB stacks exhibit distinct energy distances between the two shifted Dirac points, consistent with the polarization strength of the system. Similar results are observed during the stacking transition from the ABAB to the second saddle point (SP2) and then to the CBAC. These findings demonstrate that the doping level of the surface Gr layers can be tuned through multiple interfacial slip events. The multi-polarized states in vdW multilayers show significant potential for neuromorphic computing systems, such as neural networks⁵⁸.

In summary, our work systematically investigates interfacial polarization in vdW multilayers. Cross-layer electrostatic coupling, present in both commensurate and incommensurate multilayers, reduces interfacial spacing and enhances interfacial polarization. The overall polarization in a multilayer system results from the cumulative polarization of each interface, accounting for the piezoelectric effect. The GSFE profile for slipping an interface in a multilayer closely resembles that of a bilayer, which facilitates slip decomposition and simplifies switching analysis. The multilayer switching pathways, whether through sequential or cooperative slipping, can be evaluated by minimizing the GSFE barrier. Band structures analysis indicates that Gr doping can be effectively modulated across multiple layers, correlating with the overall interfacial polarization. These findings can apply to a range of 2D heterolayers and help rationalize the unconventional ferroelectricity observed in graphene-based moiré heterostructures.

All the calculations are performed using the density functional theory (DFT), as implemented with the Vienna ab initio Simulation Package (VASP 5.4) code^59,60. The Perdew–Burke–Ernzerhof (PBE) functional with the generalized gradient approximation (GGA) is utilized for the exchange-correlation functional⁶¹. To account for core electrons, the projector augmented wave (PAW) method is employed⁶². The DFT-D3 method with Becke-Jonson damping is used to describe the dispersion interactions^63,64. A vacuum region with a thickness of 25 Å is introduced in the models to isolate adjacent periodic images. The kinetic cutoff energy for plane-waves is set to be 500 eV. For 2D models, a 20 × 20 × 1 k-point mesh is employed for the Brillouin zone integration. All the atomic positions are fully relaxed until the force on each atom drops below 10⁻⁴eV/Å, with the energy convergence setting to 10⁻⁷eV. The polarization is evaluated using the Berry-phase method with dipole correction^45,46, divided by the surface area of the unit cell. The ferroelectric switching pathways are calculated with the climbing nudged elastic band (cNEB) method⁶⁵.