Multiscale analysis of large twist ferroelectricity and swirling dislocations in bilayer hexagonal boron nitride

Abstract

With its atomically thin structure and intrinsic ferroelectric properties, heterodeformed bilayer hexagonal boron nitride (hBN) has gained prominence in next-generation non-volatile memory applications. However, studies to date have focused almost exclusively on small-twist bilayer hBN, leaving the question of whether ferroelectricity can persist under small heterostrain and large heterodeformation entirely unexplored. In this work, we establish the crystallographic origin of ferroelectricity in bilayer hBN configurations heterodeformed relative to high-symmetry configurations such as AA-stacking and 21.786789° twisted configurations (Σ7), using Smith normal form bicrystallography. We then demonstrate out-of-plane ferroelectricity in bilayer hBN across configurations vicinal to both the AA and Σ7 stackings. Atomistic simulations reveal that AA-vicinal systems support ferroelectricity under both small twist and small strain, with polarization switching in the latter governed by the deformation of swirling dislocations rather than the straight interface dislocations seen in the former. For Σ7-vicinal systems, where existing interatomic potentials underperform particularly under extreme out-of-plane compression, we develop a density-functional-theory-informed continuum framework—the bicrystallography-informed frame-invariant multiscale (BFIM) model, which captures out-of-plane ferroelectricity in heterodeformed configurations vicinal to Σ7 stacking. Interface dislocations in these large heterodeformed bilayer configurations exhibit markedly smaller Burgers vectors compared to interface dislocations in small-twist and small-strain bilayer hBN. The BFIM model reproduces experimental results and provides a powerful, computationally efficient framework for predicting ferroelectricity in large-unit-cell heterostructures where atomistic simulations are prohibitively expensive.

---

1 Introduction

Since the discovery of ferroelectricity in Rochelle salt by Valasek,¹ ferroelectrics have emerged as an important class of materials for the development of non-volatile memory devices.^2,3 Ferroelectricity, the property of reversing the spontaneous polarization of certain materials through the application of an electric field, enables instantaneous reading/writing operations through switching polarizations.⁴ Due to their intrinsic ferroelectricity,^5,6 strong resistance to the formation of the depolarization field,⁷ and atomically thin nature, van der Waals (vdW) homo- and heterostructures—e.g., bilayer hexagonal boron nitride, bilayer molybdenum disulfide, and molybdenum disulfide-tungsten disulfide—have been recognized as promising ferroelectric materials. In particular, the tunability of ferroelectricity through controlled spatially varying relative sliding—obtained by imposing a relative twist and/or relative strain (heterostrain) between the two layers—makes vdW structures suitable for versatile applications in nano- and micro-electronic devices.^8,9 We use heterodeformation as an umbrella term to refer to relative twist, strain, or a combination of twist and strain to the bilayers.

Ferroelectricity in heterodeformed bilayer 2D materials arises from manipulating triangular domains formed during structural reconstruction of van der Waals (vdW) structures.^6,10 These structures consist of two layers with identical or distinct lattices, each containing basis atoms with different polarities, leading to unique polarizations for different vertical stackings.^5,11 Small twists and/or strains in bilayer vdW structures undergo structural reconstruction mediated by interface dislocations, forming triangular domains with alternating polarizations.^12,13‡ These domain shapes can be modified via interfacial sliding and bending of the interface dislocations under an electric field, creating a net polarization—a phenomenon known as sliding ferroelectricity.^6,9 However, current studies of sliding ferroelectricity are mostly limited to twisted bilayer 2D materials with a twist angle, θ < 2°.^9,15 Tuning of ferroelectricity through the application of heterostrain has not been extensively explored, which motivates us to explore different heterodeformations for which ferroelectricity can be observed. In this study, we consider bilayer hexagonal boron nitride (hBN) as the representative bilayer 2D material due to its analogy to graphene,¹⁶ its superior electric, chemical, and thermal properties,¹⁷ and its capability of showing sliding ferroelectricity at high operating temperatures.^6,15

Large-scale exploration of ferroelectricity across all possible heterodeformed bilayer hBN arrangements is experimentally impractical, thus motivating alternate avenues for study. Atomic-scale simulations can be very useful for such purposes, though such explorations have two limitations. First, an enormous periodic simulation domain must be designed to remove edge effects and consistently predict ferroelectricity (in the large-body limit) in a heterodeformed bilayer hBN. The required domain size can be very large even for small heterodeformations (e.g., small twist angles).¹⁸ Predicting structural reconstruction in such large domains is computationally expensive. Consequently, prior studies have largely focused on small heterodeformations, as a systematic exploration of the full heterodeformation space remains prohibitively expensive using atomistic simulations. Second, the predictability of accurate ferroelectric behavior in heterodeformed bilayer hBN requires a reliable interatomic potential to predict consistent structural reconstruction results across every heterodeformation. While structural reconstruction for small-heterodeformed bilayer hBN can be accurately predicted using atomic-scale simulations, the same cannot be achieved for large heterodeformed bilayer hBN.§ Here, we aim to develop a bicrystallography-informed frame-invariant multiscale (BFIM) model for predicting structural reconstruction in both small and large heterodeformed bilayer hBN under the presence and absence of applied out-of-plane electric fields. While previously available multiscale models are developed for the exploration of ferroelectricity in small heterodeformed bilayer 2D materials only,^8,12 the BFIM model is unique in its potential to capture ferroelectricity in large heterodeformed bilayer hBN.

The hypothesis of observation of ferroelectricity in large heterodeformed bilayer 2D materials originates from the work of Ahmed et al.¹⁹ who showed that structural reconstruction through the formation of interface dislocations is not limited to small heterodeformation in bilayer 2D materials. In this study, we show through density functional theory (DFT) simulations that two degenerate energy minima can be observed in the generalized stacking fault energy plot of 21.786789° twisted bilayer hBN, similar to 0° bilayer hBN, which is not observable through atomic-scale calculations using existing interatomic potentials. Moreover, we show that two energy minima corresponding to the 21.786789° twisted bilayer hBN demonstrate alternating out-of-plane polarization similar to the 0° bilayer hBN, which is an indicator of the sliding ferroelectricity in large heterodeformed bilayer hBN. Following this, we employ quantum-scale information of stacking energy and polarization to develop the BFIM model for efficient and effective prediction of ferroelectricity of an arbitrary heterodeformed bilayer hBN.

The remainder of this paper is organized as follows. Section 2 begins with a review of structural reconstruction and the formation of interface dislocations using molecular dynamics (MD) simulations of bilayer hBN subjected to small heterodeformations. Subsequently, MD simulations of GSFE are used to identify the translational invariances in 0° and 21.786789° twisted bilayer hBN, showing the unreliability of existing interatomic potentials to predict ferroelectricity in such configurations. In section 3, we present the BFIM model for predicting structural reconstruction in heterodeformed bilayer hBN under the presence and absence of applied electric fields. In section 4, we compare the computed results for small-heterodeformed bilayer hBN using the BFIM model with the results from atomistic simulations using LAMMPS. Following this, we show the ferroelectric domain formation in large heterodeformed bilayer hBN using the BFIM model. We summarize and conclude in section 5.

Notation: Lowercase bold letters are used to denote vectors, while uppercase bold letters represent second-order tensors unless stated otherwise. The gradient and divergence operators are denoted by the symbols ∇ and Div, respectively. We use the symbol to denote the inner product of two vectors or tensors.

2 Atomistic investigation of ferroelectricity in bilayer hBN: role of interface dislocations

Ferroelectricity refers to the reversal of spontaneous polarization in certain materials when an external electric field is applied. In bilayer hBN subjected to a small heterodeformation, ferroelectricity arises due to the presence of AB and BA domains with opposite polarizations normal to the interface. The domains are separated by interface dislocation lines,^8,9,13 which form as a result of atomic reconstruction in 2D homo/heterostructures subjected to small heterodeformation. The ferroelectric response of heterodeformed hBN results from the expansion and contraction of the domains when subjected to an out-of-plane electric field. Since the evolution of the domains is driven by the motion of interface dislocations, it is essential to understand the properties of these dislocations to quantify ferroelectricity in heterodeformed bilayer hBN.

This section aims to explore the role of interface dislocations in mediating ferroelectricity through atomistic simulations. We begin by examining the structure of interface dislocations in bilayer hBN under small relative twists and strains. Then, we demonstrate ferroelectricity by applying an electric field, which causes the dislocation to move, and leads to the expansion and contraction of the AB and BA domains. Finally, we provide arguments supporting the observation of ferroelectricity in bilayer hBN beyond the small heterodeformation range. We use Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)²⁰ to model a heterodeformed bilayer hBN at the atomic scale. The top hBN layer, represented by lattice , is constructed using the structure matrix:

where the columns of A represent the basis vectors, and a = 2.51 Å is the lattice constant of strain-free hBN. The two basis atoms are positioned at coordinates (0, 0) and relative to the basis vectors of . The bottom layer, represented by lattice , is constructed using the structure matrix B = FA, where F is the heterodeformation gradient. The bilayer is in AA stacking if F = I, and relative translations between the two layers will lead to the AB and BA stackings, as shown in Fig. 1.

Fig. 1 Three characteristics stacking configuration in bilayer hBN. The parallelepiped marked by dotted lines identifies the smallest unit cell.

To eliminate edge effects in all simulations, periodic boundary conditions (PBCs) are enforced in the plane of the bilayer. To impose PBCs, the simulation box size must be carefully selected so that the box vectors belong to the deformed lattices of the top and the bottom layers. We use the Smith normal form (SNF) bicrystallography framework developed by Admal et al.²¹ to calculate the box vectors. Additionally, SNF bicrystallography describes the translational symmetry of an interface, which will be used in section 3 to explore ferroelectricity in hBN under large heterodeformations.

The intralayer nearest-neighbor interaction between boron and nitrogen atoms (B–N bond) is modeled using the modified Tersoff potential.^22–24 The interlayer van der Waals interaction is modeled using the registry-dependent interlayer potential (ILP)²⁵ with a cutoff radius of 16 Å to ensure adequate sampling of dispersive interactions. Hexagonal boron nitride consists of boron and nitrogen atoms, which possess average partial charges of +0.4 and −0.4, respectively, making hBN a polar material.²⁶ To account for electrostatics, we included a Coulomb potential with a cutoff of 16 Å, consistent with the ILP cutoff. To simulate substrate effects in experiments, a continuum substrate is placed to interact with the hBN layers, as opposed to conventional free-suspended or out-of-plane constrained boundary conditions along the X₃ direction.^27,28 The implementation of the continuum substrate is described in section S1 of the SI. Atomic reconstruction is simulated by minimizing the total energy using the fast inertial relaxation engine (FIRE) algorithm²⁹ with an energy tolerance and force tolerance of 1 × 10⁻²⁰ eV and 1 × 10⁻²⁰ eV Å⁻¹, respectively. The resulting atomic displacements are analyzed to interpret them in terms of interface dislocations. To calculate the polarization of the system, we employ LAMMPS's built-in dipole moment calculator. Due to the non-uniqueness associated with the absolute polarization of periodic systems, we ultimately measure changes in polarization. In the next section, we discuss the characteristics of interlayer dislocations in bilayer hBN under small heterodeformations.

2.1 Characterization of interface dislocations in heterodeformed bilayer hBN

The structure of interface dislocations in hBN varies significantly depending on the relative deformation between the two layers of hBN. To demonstrate this, we conduct structural reconstruction simulations under the following two heterodeformations: (a) a 0.2992634° twist and (b) a pure stretch ofrelative to the 0° twisted AA-stacked bilayer hBN. The corresponding periodic simulation box vectors, computed from SNF bicrystallography, are shown below:

Twist:

b₁ = (480.394592e₁) Å, (2a)

b₂ = (240.197294e₁ + 416.033919e₂) Å; (2b)

Equi-biaxial strain:

b₁ = (593.619081e₁) Å, (2c)

b₂ = (293.036378e₁ + 516.248841e₂) Å, (2d)

where e_is denote the unit vectors parallel to the global axes X₁, X₂, and X₃. The PBCs ensure that the average heterodeformation stays constant during structural relaxation.

The atomic energy map of the relaxed 0.299263° twisted bilayer hBN in Fig. 2a shows the formation of two isoenergetic AB and BA domains, which are separated by interlayer dislocations. Fig. 2b shows a line scan of the displacement components along segment ①. The jumps in displacement components u₁ and u₂ in Fig. 2b suggests the Burgers vector (jump in the displacement vector) is parallel to the dislocation line. In other words, the interface dislocations in a twisted bilayer hBN have a screw character as supported by prior studies.^14,30

Fig. 2 Atomic reconstruction in a 0.299° twisted bilayer hBN. (a) Atomic energy per atom [in eV] plot showing a triangular network of interface dislocations that separate the low-energy Bernal stackings. (b) Displacement components u₁ and u₂ [in Å], measured along the line ① (in (a)) and relative to the untwisted (AA stacking) configuration, signify the screw characteristic of dislocations in a twisted bilayer hBN.

On the other hand, a 0.4219% equi-biaxial heterostrained bilayer hBN forms a spiral triangular network of dislocation lines, as illustrated in Fig. 3a. Topologically, a triangular network of straight edge dislocations is consistent with the incompatibility associated with an equi-biaxial heterostrain. However, the dislocation lines twist by 60° at the AA junctions, resulting in a swirling network. This swirling occurs because the pure edge dislocations transform locally near the AA junctions to enhance their screw character, thereby releasing the excess strain energy.^31,32

Fig. 3 Atomic reconstruction in a 0.422% equi-biaxial heterostrained bilayer hBN. (a) Atomic energy per atom [in eV] plot showing a triangular network of swirling dislocations that separate low-energy Bernal stackings. (b) Displacement components u₁ and u₂ [in Å], measured along line ① (in (a)) and relative to the undeformed (AA stacking) configuration, signify the mixed characteristics of dislocations under equi-biaxial heterostrains.

Since the formation of AB and BA domains during structural relaxation results from local relative translations between the layers, the generalized stacking fault energy (GSFE) map, which describes the interfacial energy as a function of relative translations, encodes a bilayer's dislocation properties. In particular, the periodicity of GSFE conveys the translational invariance of an interface and defines the sets of interlayer dislocations an interface can host.^19,21 Fig. 4 shows the GSFE of 0°-twist bilayer hBN calculated using density functional theory (DFT), with details in section S3 of the SI, and LAMMPS. The agreement between Fig. 4a and b confirms the accuracy of the interatomic potential used in our atomistic simulations. From the figure, we note that three other equivalent low-energy stackings surround each low-energy stacking. This feature leads to the formation of a triangular dislocation network in bilayer hBN. Additionally, the shortest distance between two degenerate minima determines the Burgers vector of interlayer dislocation, which matches the Burgers vector calculated from the plots in Fig. 2b and 3b. While the GSFE encodes the Burgers vectors of interface dislocations, it cannot describe their response to an external electric field. In the next section, we demonstrate ferroelectricity under small heterodeformations and introduce the polarization map to explain how dislocations respond to an electric field during ferroelectric transitions.

Fig. 4 Generalized stacking fault energy [meV Å⁻²] of 0° twist bilayer hBN computed using (a) DFT (b) LAMMPS, and plotted as a function of the relative displacement between the two layers. The four corners correspond to the AA stacking.

2.2 Ferroelectricity in bilayer hBN under small heterodeformations

In this section, we simulate the ferroelectric transition in a small-heterodeformed bilayer hBN and reveal its crystallographic origin via the polarization landscape. From the previous section, recall that the equi-sized AB and BA triangular domains formed during the structural relaxation of a small-heterodeformed bilayer hBN are energetically equivalent. However, their response to an applied electric field varies due to their polarity. In AB stacking, boron (positively polarized) sits above nitrogen (negatively polarized), whereas in BA stacking, the sequence is reversed. This polarity difference results in opposite reactions when an electric field is applied perpendicular to the layers.

To illustrate ferroelectricity in a small twist bilayer hBN, we subject the 0.299263° twisted hBN to an external electric field of 5 V Å⁻¹ in the ±X₃ direction. An electric field E is enforced using the fix efield command in LAMMPS, resulting in force qE on atoms with fixed charge q. The electric field induces a transformation in the bilayer structure, as shown in Fig. 5. The blue regions in the figure represent AB stacking, and the green regions indicate BA stacking. Under a positive electric field, the AB domains expand while the BA domains shrink (Fig. 5a). The bending of the interlayer dislocations mediates the transformation from one stacking to another. Next, we examine ferroelectricity in a relaxed, equi-biaxially heterostrained hBN. Recall from section 2.1 that the equi-biaxial heterostrain causes spiral dislocations. Similar to the small-twist case, the AB and BA domain areas change when an external electric field is applied, as shown in Fig. 6. However, under the out-of-plane electric field of 5 V Å⁻¹ in the −X₃ direction, the AB-stacked regions occupy ≈41% of the total simulation domain in the 0.299263° twisted bilayer hBN compared to the ≈38% in 0.4219% equi-biaxial heterostrained bilayer hBN. This corresponds to a 3% smaller AB fraction in the heterostrained bilayer, highlighting the more pronounced area change due to spiral dislocation formation.

Fig. 5 The AB (blue) and BA (green) domains in a relaxed 0.299° twisted bilayer hBN subjected to an electric field of 5 V Å⁻¹ in the (a) +X₃ direction and (b) −X₃ direction.

Fig. 6 Distribution of AB and BA regions in 0.422% equi-biaxially heterostrained bilayer hBN while subjected to applied electric field. (a) and (b) correspond to the relaxed configurations while subjected to positive and negative out-of-plane electric fields, respectively.

To predict the structural relaxation in the presence of an electric field, the GSFE alone is insufficient as the AB and BA domains are energetically equivalent. Area changes due to the electric fields are governed by the polarization landscape (PL)—defined as the polarization density as a function of relative translation between the layers. Therefore, characterizing the PL is necessary to correlate the applied electric field with the spatial variation of the stacking. Fig. 7 shows the out-of-plane polarization density map of AA-stacked bilayer hBN computed using DFT. For completeness, we also show the in-plane polarization density in section S4 of the SI, which is comparable in magnitude to the out-of-plane polarization, indicating the potential for exhibiting strong in-plane ferroelectricity in small-heterodeformed bilayer hBN.¶ For comparison, the polarization densities of AB and BA stackings relative to the AA stacking, measured using LAMMPS, are 0.264 × 10⁻³ C m⁻² and −0.265 × 10⁻³ C m⁻² respectively. Conversely, the corresponding DFT measurements from Fig. 7 are −3 × 10⁻³ C m⁻² and 3 × 10⁻³ C m⁻², which are 10 times the LAMMPS measurements with an opposite trend. We attribute this discrepancy to the assumption of localized charges in atomistic calculations of polarization, as also illustrated by Riniker et al.³³ Moreover, in our atomistic simulations, an electric field as high as 1 V Å⁻¹ is required to observe noticeable aerial changes in AB/BA sackings. On the other hand, many experimental studies^34,35 measure similar aerial changes at lower electric fields of 0.016 V Å⁻¹. Such discrepancies between atomistic simulations and experiments highlight the limitations of interatomic potentials and motivate the development of DFT-informed multiscale modeling.

Fig. 7 PL [mC m⁻²], computed using DFT, and plotted as a function of the relative displacement between the two layers of an AA-stacked bilayer hBN.

From Fig. 7, we observe that the AA stacking has a zero out-of-plane polarization, while the AB and BA stackings are oppositely polarized with equal magnitude in the out-of-plane direction. Therefore, the equi-sized AB and BA domains in a relaxed small twisted hBN with no external electric field ensure that the homostructure is unpolarized. In the presence of an out-of-plane electric field in the +X₃ direction, the AB domains (positively polarized) grow and the BA domains shrink while the AA regions stay put, which is a consequence of the work done by the external electric field, first reported by Yasuda et al.⁶ The evolution of the domains is mediated by the bending of the dislocation lines.⁹ Due to the preferential growth of the AB domains, the homostructure attains a net positive polarization.

Based on the above discussion of moiré ferroelectricity in small heterodeformed hBN, we identify the following sufficient conditions for a bilayer configuration to demonstrate ferroelectricity:

1. The heterodeformation should be small relative to a low-energy stacking up to relative translations.|| All heterodeformations introduced to this point are relative to the AA stacking.

2. The GSFE of the low-energy stacking has degenerate minima, i.e., the minima are energetically equivalent.

3. The configurations corresponding to the two minima are oppositely polarized.

The first two conditions are responsible for dislocation-mediated structural relaxation, while the last condition ensures moiré ferroelectricity. In a recent work, Ahmed et al.¹⁹ showed that an AB-stacked bilayer graphene, twisted by 21.786789°, satisfies the first two conditions, and small heterodeformations relative to this configuration result in structural reconstruction mediated by interface dislocations. In the next section, we show that the PL of the 21.786789° configuration satisfies condition 3, and proceed to explore ferroelectricity in large-twist hBN.

2.3 Ferroelectricity in large-twist bilayer hBN

In the previous section, we identified various features of GSFE and the polarization landscape of AA-stacked bilayer hBN that lead to dislocation-mediated structural relaxation and ferroelectricity in small-heterodeformed bilayer hBN. As the heterodeformation increases, structural relaxation decreases, and interpreting the relaxed structure in terms of dislocations relative to the AA stacking breaks down as the defect cores overlap.¹⁴ Interestingly, however, the authors demonstrated in an earlier work that the dislocations interpretation reemerges for small heterodeformations relative to the defect-free 21.786789° twisted bilayer hexagonal systems under out-of-plane compression. In particular, they showed that the interface energy of 21.786789° twisted bilayer graphene, minimized for local relative translations, is a local minimum with respect to the twist angle. Under small heterodeformations of this configuration, the system relaxes by nucleating interface dislocations whose Burgers vector is times the Burgers vector observed in small-heterodeformed bilayer hBN configurations. In this section, we investigate whether small heterodeformations applied to the 21.786789° twisted bilayer hBN, which itself represents a large heterodeformation relative to AA-stacked bilayer hBN, can induce ferroelectricity.

Fig. 8 shows the atomic configuration of 21.786789° twisted bilayer hBN. The supercell or the coincident site lattice (CSL) of this configuration is 7 times the unit cell of hBN and contains 28 atoms. We will refer to this large twist bilayer hBN as the Σ7 configuration. Before plotting the GSFE and PL plots, we first infer their periodicity using the SNF bicrystallography formulated by the first and last authors in Admal et al.²¹ Specifically, the GSFE and PL plots are periodic with respect to the Σ7 configuration's displacement shift complete lattice (DSCL), which is the coarsest lattice that contains lattices and . We compute the GSFE and PL of the Σ7 bilayer hBN on the domain spanned by the DSCL basis vectors b₁ = (0.948690e₁) Å and b₂ = (0.474345e₁ + 0.821590e₂) Å. Fig. 9 shows the GSFE and PL plots, computed using DFT, of the Σ7 configuration under a 28% out-of-plane compression. The two low-energy minima for the 21.786789° twisted bilayer hBN appear only under out-of-plane compression. In the absence of compression, the generalized stacking fault energy (GSFE) surface (SI S5, Fig. S2a) shows a weak energy barrier between different stackings and no neighboring degenerate minima, which are required for the formation of partial dislocations.¹⁹

Fig. 8 Atomic configuration of 21.787° twisted bilayer hBN with red and blue colors indicating boron and nitrogen atoms, respectively. (a) Top view showing a primitive unit cell of the moiré superlattice, identified by the dashed lines that consists of 28 atoms and, (b) Side view illustrating the bilayer structure and interlayer spacing.

Fig. 9 (a) GSFE [meV Å⁻²] and (b) PL [mC m⁻²] density map plots of 21.787°-twisted bilayer hBN, calculated using DFT.

The GSFE in Fig. 9a shows two degenerate minima similar to that of the 0° bilayer hBN in Fig. 4, suggesting small heterodeformations relative to the Σ7 configuration will result in interface dislocations. The size of the Burgers vectors of such dislocations is given by the distance between the minima, which is ≈0.55 Å. Moreover, analogous to the small twist case, since each minimum in Fig. 9a is surrounded by three minima, we expect the interface dislocations to form a triangular network. While in an earlier work, we demonstrated structural relaxation in heterodeformed Σ7 bilayer graphene using atomistic simulations, we are unable to repeat such a calculation for bilayer hBN in this work since the ILP potential is not applicable for large out-of-plane compressions. Fig. 10 reflects the failure of the ILP potential—not only is the magnitude of the LAMMPS-calculated GSFE inconsistent with that in Fig. 9a, but also the distribution of the extrema does not agree.

Fig. 10 The GSFE [meV Å⁻²] of a 21.787°-twisted bilayer hBN computed using LAMMPS differs significantly from its DFT counterpart in Fig. 9a.

The PL plot in Fig. 9b shows that the Σ7 stacking has a zero out-of-plane polarization density, while the stackings corresponding to the degenerate minima have opposite polarization densities of magnitude 1.35 × 10⁻⁴ C/m². Therefore, all three sufficient conditions for ferroelectricity, outlined at the end of Section 2.2, are met by large-heterodeformed bilayer hBN vicinal to the Σ7 configuration.

The potential of ferroelectricity under large heterodeformations identified in this section and the absence of reliable interatomic potentials motivate us to develop a DFT-informed multiscale framework that is computationally efficient compared to atomistic simulations and capable of capturing large-twist bilayer hBN physics. In the next section, we present the bicrystallography-informed and frame-invariant multiscale model for the prediction of ferroelectricity in arbitrarily heterodeformed bilayer hBN.

3 A multiscale model for ferroelectricity in heterodeformed hBN

In this section, we present a bicrystallography-informed frame-invariant multiscale (BFIM) model for predicting ferroelectricity in arbitrarily heterodeformed bilayer hBN. The goal of this model is to predict the structural response and the polarization density field when a bilayer hBN is subjected to an electric field. In this work, we extend the GFK model of Ahmed et al.,¹⁹ which was used to predict structural relaxation, to include polarization.

In the BFIM model, a bilayer is described as two continuum sheets. It consists of (a) a defect- and strain-free configuration called the natural configuration, which is used to construct the energy of the system, (b) a reference configuration, with respect to which displacements are measured, and (c) a deformed configuration, which represents the deformed bilayer. If a heterodeformation is vicinal to the AA stacking, then the AA stacking is chosen as the natural configuration. On the other hand, if a heterodeformation is vicinal to the Σ7 configuration, then the Σ7 configuration is chosen as the natural configuration. A small heterodeformation is introduced by uniformly deforming one of the layers, resulting in the reference configuration. The model is designed to predict the displacements from the reference configuration to the deformed configuration arising due to atomic reconstruction.

3.1 Kinematics

In the reference configuration, the two layers of the bilayer are represented by subsets Ω^ref_t and Ω^ref_b of the 2D Euclidean point space . An arbitrary material point in the bilayer is denoted by X_t or X_b, depending on whether the point belongs to the top or the bottom layer. In this study, we ignore out-of-plane displacement since ferroelectricity is primarily governed by the in-plane reconstruction.⁹ (α = t, b)** denote time-dependent deformation maps associated with the atomic reconstructions of the respective lattices. The images of the maps ϕ_t and ϕ_b constitute the deformed configuration. At t = 0, we assume ϕ_α (X_α, 0) = X_α. Therefore, ϕ_α(X_α, t)−X_α describes the displacement field associated with atomic reconstruction.

The lattice strain in the deformed configuration is measured relative to the strain-free natural configuration. If κ_α denotes the mapping between the reference and the natural configuration, and η_α denotes the mapping from the natural configuration to the deformed configuration, we have

ϕ_α = η_α°κ_α (3)

where ° denotes function composition. The inverse of κ_α is the deformation required to form the uniformly heterodeformed bilayer hBN configuration. If F_α denotes the gradient of ϕ_α, eqn (3) implies

F_α = H_αK_α where H_α: = ∇η_α and, K_α: = ∇κ_α. (4)

To construct the system's energy, we note that the total energy includes the elastic energy due to lattice strains, the interfacial van der Waals (vdW) energy, and the work done by the external electric field through its interaction with the polarization. We build each of these energy components using frame-invariant measures. A frame-invariant kinematic measure of the lattice strain is the Lagrangian strain E_α := (H^T_αH_α − I)/2. From eqn (4), E_α can be written as

E_α = (K_α^−TF^TFK_α⁻¹ − I)/2. (5)

The interfacial vdW energy and the work done by the external electric field are associated with the interface. Specifically, these energies depend on stacking, which varies spatially and depends on the relative shift between layers. Therefore, the vdW energy and polarization energy density at a point, x, in the deformed configuration are described as functions of the relative translation vector

r(x, t) = K_tX_t − K_bX_b where X_α: = ϕ_α⁻¹(x, t). (6)

Note that r is frame-invariant by design. In the next section, we will construct the constitutive law (energy and mobility) in terms of the kinematic measures C_α and r.

3.2 Constitutive law

Expressing the total energy as , we will now construct the elastic energy , van der Waals energy , and the work done by the electric field as functionals of the unknown fields ϕ_t and ϕ_b. The elastic energy is expressed in Saint Venant–Kirchhoff form aswhereis the elastic energy density of the α-th layer, and is the fourth-order isotropic elasticity tensor with lamé constants³⁶ λ = 3.5 eV Ų and μ = 7.8 eV Ų. Note that the integral in eqn (7) is over the natural configuration, with dY_α representing a unit volume in the natural configuration Ωⁿ_α. The interfacial vdW energy originates from the interaction between the layers in the region Ω_t∩Ω_b. The expression of vdW energy is constructed aswhere e_vdW is the interfacial energy density (per unit volume in the natural configuration). Note that the factor (det H_α)⁻¹ is necessary because the integration is over the deformed configuration as opposed to the natural configuration. The vdW energy density is the GSFE per unit area of the natural configuration. Therefore, it is expressed in a form that reflects the symmetry of the GSFE aswhere v_g is the strength of the GSFE, d¹ and d² are reciprocal to the two DSCL vectors that span the domain of the GSFE, and d³ = −(d¹ + d²).

The parameters v_g and c_g are obtained by comparing eqn (10) with the atomistic GSFE of AA stacking for small twist bilayer hBN, and the DFT GSFE for 28% out-of-plane compressed large-heterodeformed bilayer hBN.

The work done by the electric field originates from its interaction with the polarized interface Ω_t∩Ω_b. Therefore, the functional is constructed as

where e_pol is the work done by an external electric field E. It is of the form

e_pol(r) = χ⁻¹E·P, (12)

where P denotes the spatially varying dipole moment that depends on local stacking and is obtained from the PL plots. The out-of-plane component of the polarization density field can be represented aswhere v_f is the maximum polarization, and is the interlayer spacing between the layers. The parameter v_f is obtained by comparing eqn (13) with the atomistic PL of the AA stacking for small twist bilayer hBN, and the DFT PL for 28% out-of-plane compressed large-heterodeformed bilayer hBN.†† The conversion factor χ in eqn (12) is obtained by comparing the areal changes in AB and BA domains, predicted by the BFIM model, to those reported by Lv et al.³⁴ in 0.3° twisted bilayer hBN. Weston et al.⁹ used a similar approach to fit their model studying structural relaxation in twisted bilayer molybdenum disulfide under an applied electric field.

To simulate structural relaxation, we minimize the total energy using the gradient flow

where m ≡ 1 is the mobility associated with ϕ_αs, and δ_ϕα denotes the variation with respect to ϕ_α.

3.3 Derivation of governing equations

In this section, we derive the governing equations for the dynamics of structural relaxation in heterodeformed bilayer hBN. Computing the variational derivatives of eqn (7), (9) and (11), and substituting them into eqn (14), we obtainwhere P_α := H_α∇e_elK_α^−T, is the 2D analog of the elastic Piola–Kirchhoff stress, which measures force in Ω_α measured per unit length in Ω^ref_α. Furthermore, r_t and r_b are given by,‡‡

r_t ≈ (K_t − K_b)X_t − H_b⁻¹(ϕ_t(X_t) − ϕ_b(X_t)), (16a)

r_b ≈ (K_t − K_b)X_b − H_t⁻¹(ϕ_t(X_b) − ϕ_b(X_b)). (16b)

Details about the numerical implementation of (15) are given in section S2 of the SI. Although the BFIM model is applicable for both finite and periodic boundary conditions, we limit its implementation to PBCs in this paper because the goal is to compare bulk scale experimental observations of Lv et al.³⁴

In the next section, we will use the BFIM model to predict structural reconstruction with and without an applied electric field for any arbitrarily heterodeformed bilayer hBN.

4 Results and discussion

This section shows simulation results of the BFIM model demonstrating how structural relaxation occurs under an external electric field in various small and large heterodeformed bilayer hBNs. Additionally, we compare the structural relaxations predicted by the BFIM model to those observed in experiments under out-of-plane electric fields. We begin by first fitting χ, the only unknown parameter of the BFIM model, by simulating the structural relaxation of 0.299263° twisted bilayer hBN at E = −0.08 V Å⁻¹. Fig. 11 shows continuum simulation results of structural relaxation at E = 0, E= +0.08 V Å⁻¹, and E = −0.08 V Å⁻¹. As expected, the model predicts the triangular network of dislocations for E = 0, and domains that shrink (grow) under E = +0.08 V Å⁻¹ grow (shrink) when the sign of the electric field switches. The parameter χ is fit using a mid-point search strategy such that the areal ratio, A_AB/A_BA at E = −0.08 V Å⁻¹, matches with the experimental measurement of Lv et al.³⁴ The fitting yields χ⁻¹ = 11. Next, we validate the parametrized model by comparing (see Fig. 11d) its predictions of A_AB/A_BA under various applied electric fields with experimental measurements of Lv et al.³⁴ Fig. 11d shows good agreement between BFIM model predictions and experiment while demonstrating the non-linear dependence of A_AB/A_BA on the applied field. In contrast to Fig. 5, the AB(BA) domain increases(decreases) in area under a negative electric field, which is consistent with the PL plot of AA stacked bilayer hBN. This discrepancy in atomistic simulations is a consequence of the incorrect sign and magnitude of the LAMMPS-calculated polarization densities of AB and BA stackings, which required a significantly strong electric field to deform the interface dislocations compared to electric fields as small as 0.0475 V Å⁻¹ used in experiments. To further illustrate this field-dependent response, Fig. 12 shows the variation of the net out-of-plane dipole moment with an electric field in 0.299° twisted bilayer hBN, computed using the BFIM model. From the figure, we note that the dipole moment increases nonlinearly with the electric field with a rapid increase at lower electric field magnitudes. Such trends are consistent with those reported by Bennett et al.³⁷ in twisted bilayer molybdenum disulfide.

Fig. 11 BFIM model prediction of ferroelectric domain formation in 0.299° twisted bilayer hBN. The blue and green regions in the polarization density contour plots depict the AB and BA domains. (a), (b), and (c) correspond to E = 0, +0.08 V Å⁻¹, and −0.08 V Å⁻¹, respectively. The latter two closely resemble the atomistic simulation plots in Fig. 5. (d) Compares the areal change ratio, , illustrated in (a), with that reported by Lv et al.³⁴ under various electric fields.

Fig. 12 Variation of net out-of-plane dipole moment with applied electric field in 0.299° twisted bilayer hBN, computed using the BFIM model.

Moving on to small heterostrains, Fig. 13 shows ferroelectric domain formation in 0.4219% equi-biaxial heterostrained hBN under alternating positive and negative electric fields. From the plots, it is clear that the BFIM model accurately predicts the formation of spiral dislocations and their response to electric fields, as noted in Fig. 6.

Fig. 13 BFIM model prediction of ferroelectric domain formation in 0.422% equi-biaxial heteorostrained bilayer hBN under applied electric fields of (a) +0.08 V Å⁻¹ and (b) −0.08 V Å⁻¹.

Finally, we will use the BFIM model to demonstrate ferroelectric domain formation in large heterodeformed bilayer hBN. Recalling from section 2.3, we expect that any small heterodeformation relative to the Σ7 twisted bilayer hBN configuration will result in interface dislocation-mediated structural reconstruction and ferroelectricity. Therefore, using SNF bicrystallography, we choose 21.786789° + 0.170076° twisted bilayer hBN as a case study.

b₁ = (319.600e₁) Å, b₂ = (159.800e₁ + 276.781e₂) Å. (17)

Fig. 14a shows relaxed 21.956865° twisted bilayer hBN at 28% out-of-plane compression, consisting of a triangular network of interface dislocations. Although the dislocations appear similar to those in a small-twist bilayer hBN, the Burgers vector is markedly different. Fig. 14b shows the displacement components along the scanning direction ①. The displacement component normal to the line direction is negligible, while the component along the line direction is ≈0.55 Å, which implies that the dislocations have a screw character with Burgers vector magnitude ≈0.55 Å. Recall that this magnitude is equal to the distance between the minima of the Σ7 GSFE plotted in Fig. 9. To estimate the minimum energy required to initiate dislocation motion, we compute the Peierls stress for interface dislocations (SI S6), which is negligible, indicating that ferroelectric switching via dislocation motion occurs nearly instantaneously.

Fig. 14 (a) Relaxed configuration of 21.957° twisted bilayer hBN at 28% out-of-plane compression, computed using the BFIM model, showing an array of dislocations whose Burgers vector is shown in (b) through the displacement components along the scanning direction ① shown in (a).

From the polarization landscape in Fig. 9b, we know that the domains in Fig. 14b are oppositely polarized. Therefore, as in the small-twist case, Fig. 15 shows that the domains grow/shrink under an external electric field. However, the areal change ratio, , is ≈4 times smaller than that observed in 0.299° twisted bilayer hBN under an electric field of −0.16 V Å⁻¹. This reduction in areal change is a consequence of the reduced PL strength of large-twist bilayer hBN.

Fig. 15 BFIM model prediction of ferroelectric domain formation in 21.957° twisted bilayer hBN at 28% out-of-plane compression under applied electric fields of (a) +0.16 V Å⁻¹ and (b) −0.16 V Å⁻¹.

5 Summary and conclusions

It is well known that 2D bilayers undergo structural reconstruction when a small relative twist is applied between their layers, leading to alternating domains of low-energy stackings (AB and BA) separated by interface dislocations or strain solitons. In the case of a small-twist hBN bilayer, its two energetically equivalent AB and BA stackings are oppositely polarized. This feature results in ferroelectricity in small-twist bilayer hBN, meaning that the system becomes polarized in response to an electric field, with one domain expanding or contracting at the expense of the other, depending on the direction of the electric field.

In this paper, we studied ferroelectricity in hBN beyond the small twist case by extending previous studies to arbitrary heterodeformations. To examine the small heterostrain case, we used atomistic simulations to demonstrate ferroelectricity under small-biaxial heterostrain. In particular, we demonstrated that in this case, structural relaxation also leads to alternating AB and BA domains. However, unlike the small-twist case, they are separated by swirling interface dislocations. Alternating the direction of the electric field drives the expansion and contraction of the domains, leading to ferroelectricity.

The large-heterodeformation study was inspired by our earlier work,¹⁹ which recognized that structural reconstruction occurs not only for small heterodeformations relative to AA stacking but also for those relative to a specific twist angle of 21.786789°, which we referred to as the Σ7 configuration. Under such large heterodeformations, we demonstrated that structural reconstruction occurs through the formation of alternate low-energy stackings, separated by interface dislocations whose Burgers vector is significantly smaller than that observed in the small-twist case. This observation led us to the question—Does bilayer hBN demonstrate ferroelectricity under large heterodeformations in the vicinity of the Σ7 configuration?

The absence of a reliable interatomic potential for large heterodeformations prompted us to leapfrog the atomic scale and develop a DFT-informed continuum model—the bicrystallography-informed frame-invariant multiscale (BFIM) model—which applies to any heterodeformation. Thereafter, we systematically addressed the above question by first mapping the generalized stacking fault energy and the polarization landscape of the defect-free Σ7 configuration using DFT calculations as functions of relative translation between the layers of the Σ7 configuration. Analogous to the small-twist case, the maps demonstrated the presence of oppositely polarized low-energy stackings, indicating that heterodeformations in the vicinity of the Σ7 configuration will exhibit ferroelectricity. Using the above maps as input to the continuum model, we demonstrated ferroelectricity under large heterodeformations.

Before concluding, we acknowledge the limitations of the current study. First, the BFIM model does not account for out-of-plane displacement. Dai et al.³⁸ and Rakib et al.³⁹ noted that interlayer dislocations can transform from straight to helical lines, forming out-of-plane bulges at the AA junction. Moreover, this out-of-plane bulge might even be responsible for forming spiral dislocations in equi-biaxially heterostrained bilayer hBN, as shown by Zhang et al.⁴⁰ To incorporate out-of-plane displacement, the constitutive law of the BFIM model should include (a) a 3D GSFE⁴¹ (as opposed to the current 2D GSFE), wherein the third dimension corresponds to the interlayer spacing, and (b) bending rigidity⁴² of the constituent 2D materials. Second, our model does not predict spontaneous nonzero polarization in the absence of an applied field, nor does it predict hysteresis. While experiments⁶ report a net dipole moment under zero electric field, its origin remains unexplored. While the effect of domain-wall pinning has not been explored in the context of bilayer 2D materials, Gao et al.⁴³ demonstrated that in bulk ferroelectric PbZr_0.2Ti_0.8O₃, the local disorder created by an isolated defect can significantly alter ferroelectric switching dynamics by modifying domain-wall velocities and locally favoring one polarization state over another. These results indicate that even weak, localized energy barriers introduced by individual defects can give rise to history-dependent switching behavior. Capturing such effects in heterodeformed bilayer hBN would require the explicit inclusion of defect-induced local energy barriers in the energy landscape, which is beyond the scope of the present work and will be explored in future studies. Beyond these limitations, a natural extension of this work is to investigate in-plane ferroelectricity in heterodeformed bilayer hBN,¹³ which can be incorporated into the current model through an in-plane polarization landscape functional. We plan to explore this in future work.

Author contributions

Md Tusher Ahmed: conceptualization, methodology, software, validation, formal analysis, investigation, data curation, writing – original draft, and visualization; Chenhaoyue Wang: software, formal analysis, investigation, data curation, and writing – original draft; Amartya S. Banerjee: software, investigation, resources, writing – review & editing, supervision, and funding acquisition; Nikhil Chandra Admal: conceptualization, software, investigation, resources, writing – review & editing, supervision, project management, and funding acquisition.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting the findings of this study are provided in the supplementary information (SI). Supplementary information: implementation details of density functional theory (DFT), molecular dynamics (MD) and multiscale model, in-plane polarization density map of AA-stacked bilayer hexagonal boron nitride, GSFE plot for 21.786789° twisted bilayer hBN at equilibrium spacing, calculation of Peierls stress from the GSFE of a 21.786789° twisted bilayer hexagonal boron nitride. See DOI: https://doi.org/10.1039/d5nr04528b.

Additional datasets or scripts used for post-processing and analysis are available from the corresponding author upon reasonable request.

Acknowledgements

NCA and TA would like to acknowledge support from the National Science Foundation Grant NSF-MOMS-2239734 with S. Qidwai as the program manager. ASB and CW would like to acknowledge support through grant DE-SC0023432 funded by the U.S. Department of Energy, Office of Science, and through the UC National Laboratory Fees Research Program of the University of California, Grant Number L25CR9003. ASB and CW also acknowledge computational resource support from UCLA's Institute for Digital Research and Education (IDRE), and the National Energy Research Scientific Computing Center (NERSC awards BES-ERCAP0033206, BES-ERCAP0028072, BES-ERCAP0025168 and BES-ERCAP0025205), a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

References

1. J. Valasek, Phys. Rev., 1921, 17, 475.
2. T. Mikolajick, S. Slesazeck, M. H. Park and U. Schroeder, MRS Bull., 2018, 43, 340–346.
3. M. Zhao, G. Gou, X. Ding and J. Sun, Nanoscale, 2020, 12, 12522–12530.
4. J. Muller, T. S. Boscke, U. Schroder, S. Mueller, D. Brauhaus, U. Bottger, L. Frey and T. Mikolajick, Nano Lett., 2012, 12, 4318–4323.
5. Z. Fei, W. Zhao, T. A. Palomaki, B. Sun, M. K. Miller, Z. Zhao, J. Yan, X. Xu and D. H. Cobden, Nature, 2018, 560, 336–339.
6. K. Yasuda, X. Wang, K. Watanabe, T. Taniguchi and P. Jarillo-Herrero, Science, 2021, 372, 1458–1462.
7. R. Mehta, B. Silverman and J. Jacobs, J. Appl. Phys., 1973, 44, 3379–3385.
8. V. Enaldiev, 2D Mater., 2024, 11, 035014.
9. A. Weston, E. G. Castanon, V. Enaldiev, F. Ferreira, S. Bhattacharjee, S. Xu, H. Corte-León, Z. Wu, N. Clark and A. Summerfield, et al., Nat. Nanotechnol., 2022, 17, 390–395.
10. A. Weston, Y. Zou, V. Enaldiev, A. Summerfield, N. Clark, V. Zólyomi, A. Graham, C. Yelgel, S. Magorrian and M. Zhou, et al., Nat. Nanotechnol., 2020, 15, 592–597.
11. F. Liu, L. You, K. L. Seyler, X. Li, P. Yu, J. Lin, X. Wang, J. Zhou, H. Wang and H. He, et al., Nat. Commun., 2016, 7, 1–6.
12. V. Enaldiev, F. Ferreira, S. Magorrian and V. I. Fal'ko, 2D Mater., 2021, 8, 025030.
13. D. Bennett, G. Chaudhary, R.-J. Slager, E. Bousquet and P. Ghosez, Nat. Commun., 2023, 14, 1629.
14. E. Annevelink, H. T. Johnson and E. Ertekin, Phys. Rev. B, 2020, 102, 184107.
15. C. Woods, P. Ares, H. Nevison-Andrews, M. Holwill, R. Fabregas, F. Guinea, A. Geim, K. Novoselov, N. Walet and L. Fumagalli, Nat. Commun., 2021, 12, 347.
16. G. Yang, C. Zhu, D. Du, J. Zhu and Y. Lin, Nanoscale, 2015, 7, 14217–14231.
17. J. Yin, J. Li, Y. Hang, J. Yu, G. Tai, X. Li, Z. Zhang and W. Guo, Small, 2016, 12, 2942–2968.
18. K. Uchida, S. Furuya, J.-I. Iwata and A. Oshiyama, Phys. Rev. B, 2014, 90, 155451.
19. M. T. Ahmed, C. Wang, A. S. Banerjee and N. C. Admal, Mech. Mater., 2024, 190, 104903.
20. 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 and T. D. Nguyen, et al., Comput. Phys. Commun., 2022, 271, 108171.
21. N. C. Admal, T. Ahmed, E. Martinez and G. Po, Acta Mater., 2022, 240, 118340.
22. J. Tersoff, Phys. Rev. B, 1988, 37, 6991.
23. D. Mandelli, W. Ouyang, M. Urbakh and O. Hod, ACS Nano, 2019, 13, 7603–7609.
24. W. Ouyang, I. Azuri, D. Mandelli, A. Tkatchenko, L. Kronik, M. Urbakh and O. Hod, J. Chem. Theory Comput., 2019, 16, 666–676.
25. W. Ouyang, D. Mandelli, M. Urbakh and O. Hod, Nano Lett., 2018, 18, 6009–6016.
26. C. Y. Won and N. Aluru, J. Phys. Chem. C, 2008, 112, 1812–1818.
27. Y. Song, D. Mandelli, O. Hod, M. Urbakh, M. Ma and Q. Zheng, Nat. Mater., 2018, 17, 894–899.
28. N. P. Kazmierczak, M. Van Winkle, C. Ophus, K. C. Bustillo, S. Carr, H. G. Brown, J. Ciston, T. Taniguchi, K. Watanabe and D. K. Bediako, Nat. Mater., 2021, 20, 956–963.
29. E. Bitzek, P. Koskinen, F. Gähler, M. Moseler and P. Gumbsch, Phys. Rev. Lett., 2006, 97, 170201.
30. P. Pochet, B. C. McGuigan, J. Coraux and H. T. Johnson, Appl. Mater. Today, 2017, 9, 240–250.
31. F. Mesple, N. R. Walet, G. Trambly de Laissardière, F. Guinea, D. Došenović, H. Okuno, C. Paillet, A. Michon, C. Chapelier and V. T. Renard, Adv. Mater., 2023, 35, 2306312.
32. M. T. Ahmed, M.-k. Choi, H. T. Johnson and N. C. Admal, ACS Appl. Mater. Interfaces, 2025, 17, 30197–30211.
33. S. Riniker, J. Chem. Inf. Model., 2018, 58, 565–578.
34. M. Lv, X. Sun, Y. Chen, T. Taniguchi, K. Watanabe, M. Wu, J. Wang and J. Xue, Adv. Mater., 2022, 34, 2203990.
35. W.-C. Fan, Z. Guan, L.-Q. Wei, H.-W. Xu, W.-Y. Tong, M. Tian, N. Wan, C.-S. Yao, J.-D. Zheng and B.-B. Chen, et al., Nat. Commun., 2025, 16, 3557.
36. J. Jung, A. M. DaSilva, A. H. MacDonald and S. Adam, Nat. Commun., 2015, 6, 6308.
37. D. Bennett and B. Remez, npj 2D Mater. Appl., 2022, 6, 7.
38. S. Dai, Y. Xiang and D. J. Srolovitz, Nano Lett., 2016, 16, 5923–5927.
39. T. Rakib, P. Pochet, E. Ertekin and H. T. Johnson, Extreme Mech. Lett., 2023, 63, 102053.
40. B. Zhang, W. Qiu, X. Liao, L. He and Y. Ni, J. Mech. Phys. Solids, 2024, 105693.
41. S. Zhou, J. Han, S. Dai, J. Sun and D. J. Srolovitz, Phys. Rev. B, 2015, 92, 155438.
42. S. Dai, Y. Xiang and D. J. Srolovitz, Phys. Rev. B, 2016, 93, 085410.
43. P. Gao, C. T. Nelson, J. R. Jokisaari, S.-H. Baek, C. W. Bark, Y. Zhang, E. Wang, D. G. Schlom, C.-B. Eom and X. Pan, Nat. Commun., 2011, 2, 591.

---

Footnotes

† These authors contributed equally to this work.

‡ Disregistry in twisted bilayer graphene increases significantly as the relative twist angle is increased beyond 2.5° (ref. 14), causing overlap between neighbouring dislocations. Thus, we refer to twists greater than 2.5° as large twist heterodeformations.

§ In a previous work by Ahmed et al.,¹⁹ we showed that structural relaxation also occurs under small heterodeformations relative to the 21.786789° twisted configuration, which is defined as large heterodeformed bilayer homostructures. Interestingly, in such cases, the Burgers vector magnitude is markedly different from that observed in small heterodeformed bilayer homostructures.

¶ Our DFT-calculated polarization is in agreement with that of Bennett et al.¹³

|| We defined the energy of bilayer stacking as the minimum interface energy over all relative translations between its layers.

** t and b stand for top and bottom layers, respectively.

†† Incorporation of 28% out-of-plane compressed GSFE and PL is essential to simulate the accurate relaxed structure in the presence and absence of applied electric fields.

‡‡ See Section 4.3 in Ahmed et al.¹⁹ for fully general expressions for r_t and r_b, arguments for their approximations in (16), and boundary conditions relevant for finite systems.

---