INTERFACIAL FERROELECTRICITY BY VAN DER WAALS SLIDING

Abstract

The technology subject of the present application concerns methods and systems for manufacturing and producing stable polarized or ferroelectric layered materials.

Claims

1-45. (canceled)

46. A process for inducing polarization in a stacked multilayered diatomic hexagonal material, the process comprising orienting any two stacked layers of a diatomic hexagonal multilayered material into a stacked parallel or a nearly parallel lattice orientation to induce internal interfacial electric field normal to the layers plane at an interface between the two stacked material layers.

47. The process according to claim 46, wherein the material is hexagonal-boron-nitride (h-BN) and/or transition-metal-dichalcogenides (TMD).

48. A process for inducing ferroelectricity to a polar diatomic hexagonal multilayered material, the process comprising applying a local electric field normal to a polar diatomic hexagonal multilayered material causing domain wall sliding, to thereby induce the ferroelectricity.

49. The process according to claim 48, comprising applying a local electric field by a biased electrode or tip above the hexagonal diatomic multilayered material surface.

50. The process according to claim 48, the process comprising applying a local electric field normal to a polar crystal of a diatomic hexagonal multilayered material, to thereby cause sliding of layers in said layered material relative to each other, to provide an array of permanent and switchable polarization domains in the crystal.

51. A process for manufacturing a ferroelectric crystal, the process comprising forming or obtaining a diatomic hexagonal multilayered material having a layered stacking configuration, wherein the material layers are stacked in a parallel or a nearly parallel lattice orientation to exhibit internal interfacial electric field normal to the layer plane of the crystal and applying electric field to said layered material to induce ferroelectric properties.

52. The process according to claim 48, wherein the polar diatomic hexagonal multilayered material has one or more same or different internal interfacial polar states.

53. The process according to claim 48, for inducing ferroelectric properties to h-BN crystal or TMD crystal.

54. A polar diatomic hexagonal multilayered material having ferroelectric properties for constructing an electronic or a photoelectric or an optical device.

55. A device implementing a polar diatomic hexagonal multilayered material having ferroelectric properties.

56. The device according to claim 55 implemented in an integrated circuit.

57. The device according to claim 56, wherein the integrated circuit is an integrated circuit memory.

58. The device according to claim 55, being selected from Random Access Memory (RAM), a flash-type memory, a ferroelectric field effect transistor, a CCD multiplexer read-out system, an integrated pyroelectric detector, an integrated surface acoustic wave device, a spatial light monitor, a microwave device, a ferroelectric tunnel junction, a ferroelectric transistor, a sensor and strain sensor.

59. A device implementing a diatomic hexagonal multilayered material, the material comprising two or more material layers oriented in a parallel or a nearly parallel lattice orientation to each other, exhibiting or having internal interfacial electric field normal to the layers plane, and having pre-determined multi-switch polarization states, each of said states being determined by summing up the number of interfaces having a polarization component pointing in one direction normal to the multilayer plane, minus the number of interfaces with a polarization component pointing in the opposite direction.

60. The device according to claim 59, wherein the material is h-BN or TMD.

61. A polarization state changing device comprising at least two pre-determined multi-switch polarization states, wherein each of the at least two polarization states comprise a difference between interfaces having a polarization pointing in a first direction normal to the multilayer plane and interfaces having a polarization pointing in a second direction that is opposite the first direction; and wherein the switch includes two or more layers of selected from the group consisting of hexagonal boron nitride (h-BN) and transition-metal-dichalcogenides (TMDs), wherein the two separate pieces are configured in a stack.

62. The device according to claim 61, wherein the switch occurs by either applying an external electric field or a mechanical stimulation to shift the domain wall and locally switch the polarization.

63. The device according to claim 61, wherein the thickness of the multilayer does not affect the polarization state.

64. The device according to claim 48, wherein the device is one of a Random Access Memory (RAM), a flash-type memory, a ferroelectric field effect transistor, a CCD multiplexer read-out system, an integrated pyroelectric detector, an integrated surface acoustic wave device, a spatial light monitor, a microwave device, a ferroelectric tunnel junction, a ferroelectric transistor, a sensor, and a strain sensor.

65. The device according to claim 48, wherein the difference is determined by summing the interfaces having a polarization component pointing in the first direction and subtracting the interfaces having a polarization component pointing in the second direction.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0078] In order to better understand the subject matter that is disclosed herein and to exemplify how it may be carried out in practice, embodiments will now be described, by way of non-limiting example only, with reference to the accompanying drawings, in which:

[0079] FIGS. 1A-1D show various high symmetry interlayer stacking configurations. (A) Top view illustration of two layers. For clarity, atoms of the top layer are represented by small circles. For either group of parallel/anti-parallel twist orientations, a relative lateral shift by one lattice spacing results in a cyclic switching between three high-symmetry stacking configurations. (B) Local-registry index (LRI, defined in SM) map of two nearly parallel rigid layers with a twist angle of 0.5°. Blue regions correspond to AA stacking, whereas AB/BA stacking appear in orange (LRI=0.86). (C) Calculated LRI map after geometry relaxation of the structure presented in panel (B). Large domains of uniform untwisted AB/BA stacking appear, on the expense of the pre-optimized AA regions. The twist is accumulated in smaller AA—like regions and in the ˜10 nm wide incommensurate domain walls (bright lines), see FIG. S1 and SM for further discussion. (D) Cross-sectional illustration of two few-layered flakes (blue and light blue regions marked the top and bottom flakes, respectively) of naturally grown h-BN (AA′), which are stacked with no twist. Plus (minus) signs mark boron (nitrogen) sites. A topographical step of a single layer switches between parallel and anti-parallel stacking orientations at the interface between the two flakes. Vertical charge displacements and the resulting net polarization P_z are marked by arrows.

[0080] FIGS. 2A-2D show direct measurement of interfacial polarization. (A) Illustration of the experimental setup. An atomic force microscope is operated in Kelvin-probe mode to measure the local potential modulation, VKP, at the surface of two 3 nm thick h-BN flakes, which are stacked with a minute twist angle. (B) VKP map showing oppositely-polarized domains of AB/BA stacking (black and white), ranging in area between ˜0.01 and 1 μm2 and separated by sharp domain-walls. (C) Surface potential along the purple line marked in (B). (D) DFT calculations of the polarization P_z per unit area in two-layer (black triangles) and four-layer (blue square) finite h-BN stacks of different lateral dimensions. The out-of-plane polarization estimated based on the measured potential is shown on the right (red star).

[0081] FIG. 3 shows the dynamic flipping of polarization orientation by domain-wall sliding. Kelvin-probe maps measured consequently from left to right above a particular flake location showing domains of up (white) and down (black) polarizations. The middle image was taken after biasing the tip by a fixed DC voltage of −20 Volts and scanning it above the blue square region shown on the left-hand image. Then the tip was biased by 10 Volt and scanned again over the same region before taking the right-hand image. Consequent domain-walls positions are marked by dashed red, green, and yellow lines. Larger white (black) domains appear after positive (negative) bias scans as a result of domain-wall motions. Note that the number of domain-walls is apparently not altered.

[0082] FIGS. 4A-D show a geometric relaxation of a Moire pattern. (A) Atomic configuration of a relaxed periodic twisted bilayer h-BN with a twisted angle of θ=0.5° . The size of the entire model system is 48.9 nm×84.7 nm B and N atoms are colored by ochre and blue, respectively. Distinctly colored AB and BA domains are obtained since atoms of the upper layer hide those of the bottom layer that reside exactly below them. (B) Interlayer distance map for the relaxed twisted bilayer. (C) The energy variation during the minimization cycles applied to the model system appearing in panel (A) plotted relative to the initial energy. The first ten points represent cycles with force tolerance set to 10{circumflex over ( )}(−3) eV/Å and the last point corresponds to the final minimization step with force tolerance of 10{circumflex over ( )}(−4) eV/Å. (D) Interlayer distance (red) and local registry index (black) along the path marked by the black dashed line in panel B. The reference LRI values of the AA and AB/BA stacking modes are −1 and 0.86 (marked by the corresponding horizontal lines).

[0083] FIGS. 5A-D show topography maps of the interface of the crystals. (A, B) Surface potential and topography maps measured simultaneously on the structure presented in the main text (FIG. 2B). The top h-BN flake thickness is uniform and includes 10 layers. Topography steps in the surface of the bottom flake are marked by dashed red lines, and its total thickness at different positions is indicated. (C, D) Additional interface between a thick (>1000 layers) bottom flake and a thin (4-7 layers) top flake. Similar potential drops between the domains are observed independent of the thickness of the structures or the substrate: graphite (SiO2) in A (C) respectively.

[0084] FIGS. 6A-D show a series of Mulliken charge distribution maps. (A) Top view of the relaxed AB stacked hydrogen terminated finite bilayer h-BN flake of 1.1 nm2 contact area. Pink, blue and white spheres represent boron, nitrogen, and hydrogen atoms, respectively. (B) The Mulliken atomic charge map of (A) calculated at the B3LYP/Def2TZVP level of theory. (C) and (D) same as (A) and (B) but for the 2.9 nm2 contact area system, respectively.

[0085] FIGS. 7A-E show series of stacked bilayer and interspacing models. (A) Top view of AB (left) and AA′ (right) stacked bilayer h-BN. The top (bottom) layer atoms are marked by small (large) circles. Lattice sites that participate in the forces F_(11 (12)) in Eq. (S9) are marked by dashed green (solid red) arrows, respectively. Similar to AA′, the eclipsed atoms in AB experience both forces, however, the hollow atoms include only F_12 (twice). Also note the zero Coulomb force in the latter case due to opposite charges of yellow/blue sites. (B) Cross section of AB stacked bilayer h-BN along the dotted black line marked in (A). α(β) indicates hollow (eclipsed) sites respectively. (C) Convergence of the total interlayer Coulomb force for h=3.3 Å with the number of Bravais lattice vectors (n_m) in the summation in Eq. (S9). (D) Inter-layer spacing (h) calculated for AA′ (solid lines) and AB (dashed lines) stacked bilayer h-BN, for different fixed values of cohesion/Coulomb ratio ε/q{circumflex over ( )}2. Black lines with q=0 correspond to graphite and show smaller h for AB than AA′ stacking as expected, while orange lines with ε/q{circumflex over ( )}2=1.5 meV show the opposite, as expected for h-BN. (E) Intra-layer displacement (marked in B) as a function of ε/q{circumflex over ( )}2. The estimated value from the experiment is marked by a red star, suggesting ε/q{circumflex over ( )}2˜3 meV.

[0086] FIG. 8 shows examples of domain-wall sliding. Consequent KPFM images of the same flake location (from left to right-hand side). Between the images a biased tip (±10 V) was scanned above the region marked by a blue square. Positive tip bias resulted in domain-wall motion that increased the white domains area over the black domains and vice versa.

[0087] FIG. 9 shows the out of plane polarization of an AB stacked multi-layer system. The out-of-plane polarization as a function of number of layers for the fully AB stacked multilayer system (squares) compared to the results presented in FIG. 10A (circles). The polarization was calculated using Gaussian.

[0088] FIGS. 10A-10C provide polarization in TMDs. (A) Potential map at the surface of three WSe2 layers stacked in a parallel orientation. Five fixed potential values are measured corresponding to the five different orientation configurations of two independent polarizations. Scale bar is 1 μm. (B) The potential difference between oppositely polarized domains is 120 mV and 240 mV for one and two interfaces respectively (dashed line cuts as shown in A). (C) DFT calculation of the potential along the cross section the two interfaces. The measured and calculated values for single and double interfaces confirm the linear polarization increase with the number of layers.

DETAILED DESCRIPTION OF THE INVENTION

[0089] As disclosed herein, the invention provides a process for inducing polarization in a stacked multilayered diatomic hexagonal material, the process comprising orienting any two stacked layers of a diatomic hexagonal multilayered material into a stacked parallel lattice orientation to induce internal interfacial electric field normal to the layers plane at an interface between the two stacked material layers.

[0090] In some embodiments, the process comprises layering two or more flakes of a diatomic hexagonal material such that a stacked multilayered material is obtained wherein each layer in the multilayered material is in a parallel lattice orientation.

[0091] In some embodiments, the process comprises depositing a layer of a diatomic hexagonal material on top a layer of same diatomic hexagonal material such that a stacked multilayered material is obtained wherein each layer in the multilayered material is in a parallel lattice orientation.

[0092] In some embodiments, the process comprises providing a diatomic hexagonal multilayered crystal or material, separating layers making up the crystal or material into separate layers and stacking at least one of the separated layers on top of another, such that an internal interfacial electric field normal to the layers plane of the material is formed at an interface between the two stacked material layers.

[0093] In some embodiments, the process comprises [0094] providing two or more flakes of a diatomic hexagonal material; [0095] stacking the two or more flakes into a multilayered material, wherein each layer in the multilayered material is in a parallel lattice orientation.

[0096] In some embodiments, the two or more flakes are obtained from exfoliating layers of a diatomic hexagonal multilayered crystal or material.

[0097] In some embodiments, the number of layers in the multilayered material is at least two.

[0098] In some embodiments, the number of layers in the multilayered material is two or three or more.

[0099] In some embodiments, the number of polarization states is equal to the number of internal interfaces between stacked layers in the multilayered material.

[0100] In some embodiments, the number of polarization states is one or is at least two.

[0101] In some embodiments, the process is for obtaining a polarized diatomic hexagonal multilayered material having one or two or more same or different internal interfacial polarized states.

[0102] In some embodiments, the diatomic hexagonal multilayered material is selected from hexagonal-boron-nitride (h-BN), transition-metal-dichalcogenides (TMD), hexagonal-aluminum-nitride (h-AlN), hexagonal-zinc-oxide (h- ZnO) and hexagonal-gallium-nitride (h-GaN).

[0103] In some embodiments, the material is hexagonal-boron-nitride (h-BN) or transition-metal-dichalcogenides (TMD).

[0104] In some embodiments, the TMD is selected from MoS.sub.2, WS.sub.2, MoSe.sub.2 and WSe.sub.2.

[0105] In some embodiments, the process is for inducing polarization in a h-BN crystal.

[0106] In some embodiments, the process is for inducing polarization in a TMD crystal.

[0107] Also provided is a diatomic hexagonal multilayered material, the material comprising two or more material layers oriented in a parallel lattice orientation to each other, exhibiting or having internal interfacial electric field normal to the layers plane.

[0108] In some embodiments, the material has one or more internal interfacial electric fields (or polarization domains), each being at a different interface between two stacked layers.

[0109] Further provided is a multipolarized diatomic hexagonal multilayered material, the material comprising two or more material layers oriented in a parallel lattice orientation to each other and two or more polarization domains.

[0110] In some embodiments, the diatomic hexagonal multilayered material is selected from hexagonal-boron-nitride (h-BN), transition-metal-dichalcogenides (TMD), hexagonal-aluminum-nitride (h- AlN), hexagonal-zinc-oxide (h- ZnO) and hexagonal-gallium-nitride (h- GaN).

[0111] In some embodiments, the material is hexagonal-boron-nitride (h-BN) or transition-metal-dichalcogenides (TMD).

[0112] In some embodiments, the TMD is selected from MoS.sub.2, WS.sub.2, MoSe.sub.2 and WSe.sub.2.

[0113] In some embodiments, the material is stable at room temperature.

[0114] In some embodiments, the material is prepared by a process according to processes of the invention.

[0115] In some embodiments, the material has three or more material layers, a magnitude of polarization of the material is equal to a polarization measured or determined for a single interface between two stacked layers multiplied by the number of layer interfaces.

[0116] In some embodiments, the material is characterized by a linear enhancement of polarization, wherein the magnitude of polarization being dependent on the number of interfaces between any two stacked layers.

[0117] In some embodiments, the material is characterized by diminished depolarization surface effects.

[0118] The invention further provides a process for inducing ferroelectricity to a polarized diatomic hexagonal multilayered material, the process comprising applying a local electric field normal to a polarized diatomic hexagonal multilayered material causing domain wall sliding, to thereby induce the ferroelectricity.

[0119] In some embodiments, the process comprises applying a local electric field by an electrode or a biased tip above the hexagonal diatomic multilayered material surface.

[0120] In some embodiments, the process is for manufacturing a room temperature stable ferroelectric crystal, the process comprising applying a local electric field normal to a polarized crystal of a diatomic hexagonal multilayered material, to thereby cause sliding of layers in said layered material relative to each other, to provide an array of permanent and switchable polarization domains in the crystal.

[0121] A process is also provided for manufacturing a room temperature stable ferroelectric crystal, the process comprising forming or obtaining a diatomic hexagonal multilayered material having a layered stacking configuration, wherein the material layers are stacked in a parallel lattice orientation to exhibit internal interfacial electric field normal to the layer plane of the crystal and applying electric field to said layered material to induce room temperature stable ferroelectric properties.

[0122] In some embodiments, the number of layers in the multilayered material is at least two.

[0123] In some embodiments, the number of layers in the multilayered material is two or three or more.

[0124] In some embodiments, the polarized diatomic hexagonal multilayered material has one or more same or different internal interfacial polarized states.

[0125] In some embodiments, the process is for inducing ferroelectric properties to h-BN crystal.

[0126] In some embodiments, the process is for inducing ferroelectric properties to TMD crystal.

[0127] Also provided is a diatomic hexagonal multilayered material having room temperature stable ferroelectric properties.

[0128] In some embodiments, the material is h-BN or TMD.

[0129] In some embodiments, the material is prepared by a process according to the invention.

[0130] A use is provided of a ferroelectric material according to the invention in constructing an electronic or a photoelectric or an optical device.

[0131] A device may be any such device implementing a ferroelectric material according to the invention. The device may be selected from integrated circuits, such as an integrated circuit memory.

[0132] In some embodiments, the device is selected from Dynamic Random Access Memory (DRAM), a flash-type memory, a ferroelectric field effect transistor, a CCD multiplexer read-out system, an integrated pyroelectric detector, an integrated surface acoustic wave device, a spatial light monitor, a microwave device, a ferroelectric tunnel junction, a ferroelectric transistor, a sensor and strain sensor.

[0133] A device implementing a material according to the invention may have pre-determined multi-switch polarization states, each of said states being determined by summing up the number of interfaces having a polarization pointing in one direction normal to the multilayer plane, minus the number of interfaces with a polarization pointing in the opposite direction.

[0134] In the following description, numerous details are set forth for the purpose of explanation. However, one of ordinary skill in the art will realize that the invention may be practiced without the use of these specific details.

Materials and Methods

a) Device Fabrication

[0135] h-BN flakes of various thicknesses (1-5 nm) were exfoliated onto a SiO.sub.2 surface. A particular flake is selected to have several topographic steps of a few-layers thickness. The flake is ripped off into two pieces which are stacked together by a polymer stamp. During the stamping processes we make sure to minimize any twist orientation. Finally, the two-flakes-sandwich is placed on a conducting graphite flake or alternatively on a gold substrate, using the same dry transfer method.

b) AFM Scans

[0136] Topography and Kelvin probe force microscopy (KPFM) measurements are acquired simultaneously (FIG. 5), using Park System NX10 AFM in a non-contact scanning mode. The electrostatic signal is measured in the first harmonic using a built-in lock in amplifier. We use Multi-75g and PPP-EFM n-doped tips with conductive coating. The mechanical resonance frequency of the tips is 75 kHz, and the force constant is 3 N/m. The cantilever is oscillated mechanically with an amplitude of ˜20 nm. The cantilever is also excited with an AC voltage to perform KPFM measurements as described below, with amplitude of 3-6 V and frequency of 17 kHz. The DC voltage is controlled by a servo motor to obtain the surface potential measurements. The images are acquired using Park SmartScan software and the data analysis is performed with Gwyddion program.

[0137] To switch the domain orientation biased scans are performed in a pin-point mode. Here the tip approaches the surface vertically at each pixel in the scanned area. The estimated maximum force during this approach is 50 nN. This mode minimizes lateral forces between the tip and the surface.

c) Kelvin-Probe Surface Potential Measurements

[0138] The AFM tip and the sample are treated as a parallel-plate capacitor model. The charge induced on the tip and the substrate is affected by the voltage applied between them, and potential drops related to the sample.

[0139] The applied voltage on the tip consists of DC and AC components. The total voltage is given by:

V=V.sub.DC+V.sub.AC sin(ωt)+V.sub.CPD

where V.sub.CPD is the contact potential difference, which originates from the different work function of the tip and the substrate. The force acting on the tip is:

[00002]$F=\frac{A}{2ϵ}{\rho }^2,$

where A is the effective area of the capacitor, ϵ is the dielectric constant and ρ is the two-dimensional charge density. The latter can be extracted from:

[00003]$V=\frac{\rho }{ϵ}d+V_{\mathrm{int}},$

where d is the distance between the plates, and V.sub.int is the voltage drop at the h-BN interface. This claim holds assuming the sample is neutral and the field outside the sample from the charges distribution in the sample is zero. After inserting it in the force equation, we get the first harmonic of the force:

F(ω)∝(V.sub.int−V.sub.DC−V.sub.CPD).sub.AC sin(ωt)

[0140] It vanishes for V.sub.DC=V.sub.int−V.sub.CPD. The main principle of KPFM is to apply a DC voltage that nullifies the first harmonic, so V.sub.int signal can be extracted from variation in the KPFM signal, V.sub.DC, at different lateral positions above the surface.

d) Model System and Classical Force-Field Calculations

[0141] To study the structural properties of twisted h-BN interfaces we constructed a model system consisting of two h-BN layers with an interlayer misfit angle of ˜0.5°. To mimic the experimental scenario, a laterally periodic supercell was constructed with a triangular lattice of periodicity L=|n{right arrow over (a)}.sub.1+m{right arrow over (a)}.sub.2|, where the primitive lattice vectors are given by {right arrow over (a)}.sub.1=a.sub.hBN(√{square root over (3,)} 0) and

[00004]${\stackrel{.fwdarw.}{a}}_2=\frac{a_{hBN}}{2}\left(\sqrt{3},3\right)$

a.sub.hBN2.505 Å based on the Tersoff potential equilibrium bond-length of b.sub.BN=1.446 Å. The indices n=195 and m=1 were chosen to fulfil the condition:

[00005]$\begin{array}{cc}\cos \left(\theta \right)=\frac{2n^2-m^2+2\mathrm{nm}}{2\left(n^2+m^2+\mathrm{nm}\right)}.& \left(\mathrm{S1}\right)\end{array}$

[0142] The corresponding moiré pattern dimension is

[00006]$L=\frac{b_{BN}}{\sqrt{2-2\cos \left(\theta \right)}}=16.3\mathrm{nm}.$

The parallelepiped supercell was then multiplicated to construct a rectangular supercell consisting of more than 300,000 atoms.

[0143] The structural properties of the twisted h-BN interface were calculated using the Tersoff intra-layer potential in conjunction with the recently developed dedicated interlayer potential (ILP). We first optimized the geometry of the top layer atoms with fixed supercell size using the Fire algorithm and keeping the bottom layer rigid. This was followed by optimization of the supercell dimensions by the conjugate gradient (CG) algorithm while scaling the rigid bottom layer according to the simulation box size. This two-step energy minimization procedure was repeated for ten times, which is sufficient to obtain well converged results (see FIG. 4C). In each repetition, both minimization stages were terminated when the forces acting on each degree of freedom reduced below 10.sup.−3 eV/Å. The system was further relaxed by Fire algorithm at a force tolerance of 10.sup.−4 eV/Å. The optimized structure exhibits atomic reconstruction with distinct AB and BA stacking domains separated by domain walls (see FIG. 4A).

e) Local Registry Index Analysis

[0144] The local registry index (LRI) (FIG. 1D) is a method introduced to quantify the degree of local interfacial registry matching at rigid material interfaces. The idea is to assign a number between −1 and 1 to each atom in the layer, signifying whether it resides in an optimal or worse stacking region, respectively. To this end, a circle is associated with each atomic position in the two layers and the overlaps between circles of one layer and those of the adjacent layer are evaluated. For the case of h-BN three types of overlaps are considered, namely S.sub.i.sup.NN, S.sub.i.sup.NB.sub.i=S.sub.i.sup.BN, and S.sup.BB.sub.i. Here, S.sub.i.sup.JKsignifies the overlap of the circle associated with atom i of type J in one layer with all circles associated with K type atoms in the adjacent layer. The radius of the circles associated with B and N atoms is taken as r.sub.B=0.15b.sub.hBN, and r.sub.N=0.5b.sub.nBN, which provides good qualitative agreement between registry index maps and the sliding potential energy surfaces obtained from density functional theory calculations. The LRI of atom i is then defined as the average registry index of itself and its three nearest neighbors (j, k, l) within the entire layer, as follows:

[00007]$\begin{array}{cc}{\mathrm{LRI}}_i=\frac{1}{3}{.Math.}_{n=j,k,l}\frac{\left[\left(S_i^{NN}+S_n^{NN}\right)-\left(S_i^{\mathrm{NN},\mathrm{opt}}+S_n^{\mathrm{NN},\mathrm{opt}}\right)\right]+\left[\left(S_i^{BB}+S_n^{BB}\right)-\left(S_i^{\mathrm{BB},\mathrm{opt}}+S_n^{\mathrm{BB},\mathrm{opt}}\right)\right]-\left[\left(S_i^{NB}+S_n^{NB}\right)-\left(S_i^{\mathrm{NB},\mathrm{opt}}+S_n^{\mathrm{NB},\mathrm{opt}}\right)\right]}{\left[\left(S_i^{\mathrm{NN},\mathrm{worst}}+S_n^{\mathrm{NN},\mathrm{worst}}\right)-\left(S_i^{\mathrm{NN},\mathrm{opt}}+S_n^{\mathrm{NN},\mathrm{opt}}\right)\right]+\left[\left(S_i^{\mathrm{NN},\mathrm{worst}}+S_n^{\mathrm{NN},\mathrm{worst}}\right)-\left(S_i^{\mathrm{NN},\mathrm{opt}}+S_n^{\mathrm{NN},\mathrm{opt}}\right)\right]+\left[\left(S_i^{\mathrm{NB},\mathrm{worst}}+S_n^{\mathrm{NB},\mathrm{worst}}\right)-\left(S_i^{\mathrm{NB},\mathrm{opt}}+S_n^{\mathrm{NB},\mathrm{opt}}\right)\right]}& \left(\mathrm{S2}\right)\end{array}$

where S.sub.i.sup.JK,opt and S.sub.i.sup.JK,worst are S.sub.i.sup.JK evaluated at the optimal and worst local stacking modes, respectively (AA′ and AA in the case of h-BN, respectively, see FIG. 1A). The calculated LRI.sub.i is then transformed by −(2LR.sub.i−1) to make it range be between [−1, 1]. With this, the LRI at an AA′ (AA) stacked region is 1 (−1) respectively, and that of an AB stacked region is 0.86.

[0145] Plotting the LRI following geometry relaxation as discussed above (FIG. 1D) we found an ordered array of AB and BA stacked domains separated by sharp domain walls. To estimate the width of the domain-wall region, we plotted a cross section of the out-of-plane height profile (see FIG. 4D) along the path marked by the dashed black line in FIG. 5B. A clear inverse correlation between the height map (red line) and the registry index (black line) is obtained. At the AB and BA stacked regions the interlayer distance is relatively constant at ˜3.23 Å and correspondingly a relative constant value of LR˜0.86 is obtained. At the center of the domain wall, the interlayer distance increases by ˜0.05 Å and the LRI reduces to ˜0.67, whereas at the domain wall crossings the interlayer distance reaches ˜3.57 Å and the LRI drops to ˜−0.93. Using a Gaussian fitting to the height profile near the domain wall we can estimate the domain wall width to be ˜10 nm.

f) Dipole Moment Calculations

[0146] To evaluate the dipole moment developing in the system we considered a finite AB stacked hexagonal h-BN bilayer model with a surface area of 1.1 nm.sup.2 and armchair edges. The flake was initially constructed with uniform B—N bond lengths of 1.446 Å and the edges were saturated by hydrogen atoms with initial B—H and N—H bond lengths of 1.200 Å and 1.020 Å, respectively (FIG. 6A,C). The structures were optimized using the hybrid B3LYP exchange-correlation density functional approximation and the double-ζ polarized 6-31G** Gaussian basis set as implemented in the Gaussian 16 suite of programs. This was followed by refined relaxation adding Grimme's D3 dispersion correction and using the 6-31+G** basis set. Finally, single point calculations were performed on the minimized structures at the B3LYP/6-31+G** and B3LYP/Def2TZVP level of theory. Comparison of the out-of-plane dipole moment components obtained using the three basis sets is provided in Table 1, showing that our results are well converged with respect to basis set size. The value calculated by Def2TZVP was used in FIG. 1C in the main text.

[0147] To verify that the flake size used is sufficiently large, we repeated the dipole moment calculation for a bilayer flake with surface area of 2.9 nm.sup.2. As seen in Table 1 the obtained values are within 20-25% with those of the smaller flake, indicating that edge effects are relatively small and validating the qualitative value of the results. This is further supported by the Mulliken charge analysis map provided in FIG. 6B, D, showing that most of the charge is located in the bulk region of the flakes.

TABLE-US-00001 TABLE 1 The out-of-plane dipole moments (p.sub.z) obtained for two AB stacked bilayer h-BN finite flake models. Dipole moment (Debye/nm.sup.2) Area (nm.sup.2) 6-31G** 6-31 + G** Def2TZVP 1.1 0.66 0.55 0.55 2.9 0.52 0.45 0.45

g) Minimalistic Classical Cohesive Energy Model

[0148] Our minimalistic model provides a classical estimate of the out-of-plane polarization of the AB bilayer interface treating the boron and nitrogen atoms as point charges (see FIG. S4A,B), interacting via Pauli and van der Waals (VdW) forces (described by the Lennard-Jones (LJ) potential), and Coulomb interactions. The total interlayer energy is written as follows:

[00008]$\begin{array}{cc}E=\frac{1}{2}{.Math.}_{i,j}\left[4\epsilon \left({\left(\frac{\sigma }{r_{ij}}\right)}^{12}-{\left(\frac{\sigma }{r_{ij}}\right)}^6\right)+\frac{q_i{q}_j}{r_{ij}}\right],& \left(\mathrm{S3}\right)\end{array}$

[0149] where ε is the cohesive energy with σ≡3.3 Å. We note that realistic models of h-BN should take atom specific ε and σ values. Here, however we are interested in a qualitative description of the system and hence, for simplicity, we limit the treatment for uniform parameter values. The differences in electronegativity of the boron and nitrogen atoms are effectively taken into account by assigning dimensionless partial charges located at the nuclear centers q=±q.sub.i/e for i∈B, N respectively. The parameter q.sup.2/εσ controls the relative strength between the Coulomb and LJ interactions. As we will demonstrate, this competition determines the sign of the polarization at the AB interface. We denote by α the atomic sites in one layer that reside above hexagon centers in the other layer (termed herein as hollow sites). Correspondingly, β denotes atomic sites in one layer that reside above oppositely charged sites on the adjacent layer (termed herein as eclipsed sites). In each layer we use h.sub.α or h.sub.β to denote vertical heights of α and β atomic sites, measured with respect to the midplane of the AB interface. To compute the polarization, we minimize the classical energy with respect to h.sub.α, h.sub.62 , via an approximate two-step protocol: [0150] 1. First, we set h.sub.α=h.sub.62 =h/2 and minimize the interaction energy with respect to h. [0151] 2. Then we allow for finite relative vertical motion of B-N pairs around the optimal h value 2Δd=h.sub.α−h.sub.β, which generates the polarization. Note that no lateral atomic motion is allowed.

Optimal Interlayer Spacing at the AA′ Stacking Mode

[0152] As a reference, we first consider two h-BN layers in the AA′ stacking configurations with h.sub.α=h.sub.β=h/2. The total force per atom is:

[00009]$\begin{array}{cc}{F}_{\mathrm{AA}\text{'}}\left(h\right)=-\frac{{\mathrm{dE}}_{{\mathrm{AA}}^{\prime }}\left(h\right)}{\mathrm{dh}}=F_{\mathrm{AA}\text{'}}^{\mathrm{LJ}}+F_{\mathrm{AA}\text{'}}^C& \left(\mathrm{S4}\right)\end{array}$

[0153] The Coulomb contribution can be written as F.sub.AA′.sup.C=−F.sub.11.sup.C+F.sub.12.sup.C (with F.sub.11.sup.C, F.sub.12.sup.C>0), where

[00010]$\begin{array}{cc}\begin{array}{cc}{F}_{11}^C\left(h\right)={.Math.}_{{\stackrel{.fwdarw.}{R}}_{11}}\frac{e^{2}h}{{\left({\stackrel{.fwdarw.}{R}}_{11}^2+h^2\right)}^{3/2}},& F_{12}^C\left(h\right)={.Math.}_{{\stackrel{.fwdarw.}{R}}_{12}}\frac{e^{2}h}{{\left({\stackrel{.fwdarw.}{R}}_{12}^2+h^2\right)}^{3/2}}\end{array}.& \left(\mathrm{S5}\right)\end{array}$

[0154] Here, {right arrow over (R)}.sub.11={right arrow over (R)}.sub.n1,n2 denote in plane lattice vectors connecting equivalent atoms, namely Bravais lattice vectors, and {right arrow over (R)}.sub.12={right arrow over (R)}.sub.n1,n2−{circumflex over (x)}R.sub.0 denote in plane lattice vectors connecting inequivalent atoms, where {circumflex over (x)}R.sub.0 is a vector connecting nearest-neighbors. The corresponding Bravais lattice vectors of the honeycomb lattice are given by {right arrow over (R)}.sub.n.sub.1.sub.,n.sub.2=n.sub.1{right arrow over (R)}.sub.1+n.sub.2{right arrow over (R)}.sub.2, and

[00011]${\stackrel{.fwdarw.}{R}}_{1,2}=R_0\left(\frac{3}{2},\pm \frac{\sqrt{3}}{2}\right),$

with R.sub.0=1.4 Å. Quick convergence of F.sub.AA′.sup.C, is guaranteed if for any pair of integers n.sub.1, n.sub.2 the term F.sub.11.sup.C(h) is combined with F.sub.12.sup.C(h) calculated for −n.sub.1, −n.sub.2 and the sums are taken over the range −n.sub.m≤n.sub.1, n.sub.2≤n.sub.m with sufficiently large n.sub.m. The force then converges as 1/n.sub.m.sup.2 (not shown) in FIG. S4C. Notably, the attractive force e.sup.2/h.sup.2 associated with a single vertical bond is strongly suppressed due to the alternating charges within the layer and the small ratio R.sub.0/h. This reduces the bare Coulomb interlayer energy

[00012]$\frac{e^2}{h}\sim 4.3\mathrm{eV}$

into the meV regime, comparable with the VdW scale ε.sup.43.

[0155] Similarly, the LJ force can be split as F.sub.AA′.sup.LJ=F.sub.11.sup.LJ+F.sub.12.sup.LJ, where

[00013]$\begin{array}{cc}{F}_{11}^{LJ}\left(h\right)=4\epsilon {.Math.}_{{\stackrel{.fwdarw.}{R}}_{11}}\left(\frac{12{\sigma }^{12}h}{{\left({\stackrel{.fwdarw.}{R}}_{11}^2+h^2\right)}^7}-\frac{6{\sigma }^{6}h}{{\left({\stackrel{.fwdarw.}{R}}_{11}^2+h^2\right)}^4}\right),F_{12}^{LJ}\left(h\right)=4\epsilon {.Math.}_{{\stackrel{.fwdarw.}{R}}_{12}}\left(\frac{12{\sigma }^{12}h}{{\left({\stackrel{.fwdarw.}{R}}_{12}^2+h^2\right)}^7}-\frac{6{\sigma }^{6}h}{{\left({\stackrel{.fwdarw.}{R}}_{12}^2+h^2\right)}^4}\right).& \left(\mathrm{S6}\right)\end{array}$

[0156] The zero-force condition yields the optimal interlayer distance h, marked in FIG. 7D by solid lines. As shown h decreases upon decreasing the Lenard-Jones energy scale ε with respect to the Coulomb energy.

[0157] We note that a reasonable approximation for F.sub.AA′.sup.LJ(h), for small R.sub.0/h, consists of treating the particles as a uniform mass distribution, i.e. replacing the sum over {right arrow over (R)} by integration, yielding

[00014]$\begin{array}{cc}{F}_{\mathrm{continuu𝔪}}^{\mathrm{LJ}}\left(h\right)=\underset{R_0/h.fwdarw.0}{\lim }{F}^{\mathrm{LJ}}=\epsilon \frac{32\pi }{3\sqrt{3}{R}_0^2}h\left(\frac{{\sigma }^{12}}{h^{12}}-\frac{{\sigma }^6}{h^6}\right).& \left(\mathrm{S7}\right)\end{array}$

Optimal Interlayer Spacing at the AB Stacking Mode

[0158] Now consider two h-BN layers at the AB stacking configuration. The corresponding interlayer force can be written as F.sub.AB=F.sub.AB.sup.LJ+F.sub.AB.sup.C. The unit cell consists of two types of atomic sites, one type where atoms of the two layers reside atop of each other (eclipsed) and the other type where an atom of one layer resides atop a hexagon center of the other layer (hollow sites). Note that the Coulomb contribution of the hollow sites vanishes due to symmetry considerations. Hence, we only have the Coulomb contribution from the eclipsed atomic sites. Since the latter have the exact same configuration in the AA′ and AB stacking modes (see FIG. 7A) the overall Coulomb force contribution per atom in the AB stacking mode is half of that in the AA′ mode,

[00015]$F_{AB}^C=\frac{1}{2}{F}_{A{A}^{\prime }}^C.$

[0159] For the VdW part, similar to the AA′ stacking case, the eclipsed atomic sites give F.sub.11.sup.LJ+F.sub.12.sup.LJ of Eq. S10, whereas the hollow atomic sites give 2F.sub.12.sup.LJ due to the unique symmetry of the AB stacked bilayer hexagonal lattice. Therefore, in total we obtain for the LJ force contribution per atom that

[00016]$F_{AB}^{LJ}=\frac{F_{11}^{LJ}+F_{12}^{LJ}}{2}+\frac{2F_{12}^{LJ}}{2}=\frac{F_{11}^{LJ}}{2}+\frac{3F_{12}^{LJ}}{2}.$

The zero-force condition yields the optimal AB stacking interlayer distance, marked by the dashed gray line in FIG. 7D. Similar to the case of AA′ stacking mode, h decreases upon decreasing E/q.sup.2. Note however, that while h (AA′)>h (AB) for large ε/q.sup.2 the situation is inverted for

[00017]$\frac{\epsilon }{q^2}\lesssim$

1.5 meV. Specifically, when q.fwdarw.0 our model corresponds to the case of graphite with optimal AB stacking mode.

Relative Vertical Displacement

[0160] On top of the interlayer spacing, we now allow a small opposite motion of the hollow site (α) and eclipsed site (β) atoms: h.sub.a=h/2+Δd,h.sub.β=h/2−Δd, (S8) as marked in FIG. S4B. We now analyze the total energy of the system as a function of Δd, via the Harmonic approximation

[00018]$E\cong \mathrm{const}+{\left(\frac{\mathrm{dE}}{d\left(\Delta d\right)}\right)}_{\Delta d=0}\Delta d+\frac{1}{2}{\left(\frac{d^{2}E}{{d\left(\Delta d\right)}^2}\right)}_{\Delta d=0}$

Δd.sup.2 around Δd.sup.2=0. The assumption that Δd<<R.sub.0 can be justified by the experimental estimate of Δd˜10.sup.−3 Å (see main text). Our simple model considered here, similar to an Einstein model for lattice vibrations, assumes that the Δd coordinates are independent Harmonic oscillators with spring constant

[00019]$K_{\Delta d}=\frac{d^{2}E}{{d\left(\Delta d\right)}^2}.$

Crucially, Δd=0 is not a minimum due to the reduced symmetry of the AB interface, with alternating eclipsed and hollow sites, which imposes a finite normal force (per atom) of:

[00020]$\begin{array}{cc}{F}_{\Delta d}=-\frac{\mathrm{dE}}{d\left(\Delta d\right)}=F_{\Delta d}^C+F_{\Delta d}^{\mathrm{LJ}},& \left(S9\right)\end{array}$

where F.sub.Δd.sup.C and F.sub.Δd.sup.LJ are defined as the relative displacement force contributions of the Coulomb and LJ terms. Note that F.sub.Δd is the difference between the forces acting on the eclipsed versus the hollow sites.

[0161] Since, for Δd=0, hollow site atoms see a locally charge-neutral configuration on the other layer the Coulomb part of this linear force originates only from the eclipsed atomic sites. For the latter, the atom in the upper layer is attracted to the one residing exactly below it in the other layer. However, it is repelled by its next nearest neighbors in the other layer, and so on. When performing the entire lattice sum for the eclipsed site we obtain that the overall force is always attractive:

[00021]$F_{\Delta d}^C=\frac{F_{11}^C\left(h\right)-F_{12}^C\left(h\right)}{2}>0.$

(S10)

[0162] On the other hand, the LJ potential, which can be written as the difference between the hollow site (2F.sub.12.sup.LJ(h)) and the eclipsed site (F.sub.11.sup.LJ(h)+F.sub.12.sup.LJ(h)) contributions yields a repulsive force per atom near the equilibrium interlayer distance:

[00022]$\begin{array}{cc}{F}_{\Delta d}^{LJ}=\frac{2F_{12}^{LJ}\left(h\right)}{2}-\frac{F_{11}^{LJ}\left(h\right)+F_{12}^{LJ}\left(h\right)}{2}=\frac{F_{12}^{LJ}\left(h\right)-F_{11}^{LJ}\left(h\right)}{2}<0.& \left(S11\right)\end{array}$

[0163] This signifies that at the equilibrium interlayer distance, the eclipsed site contribution is more repulsive than that of the hollow site counterpart, mainly due to the fact that the eclipsed atoms are forced to reside within the steep Pauli repulsion wall side of their pairwise interaction.

[0164] Overall, from Eqs. S15 and S16 our crude estimate yields

[00023]$\begin{array}{cc}\Delta d\sim \frac{F_{\Delta d}}{K_{\Delta d}}=\frac{\left(F_{11}^C\left(h\right)-F_{12}^C\left(h\right)\right)+\left(F_{12}^{LJ}\left(h\right)-F_{11}^{LJ}\left(h\right)\right)}{2K_{\Delta d}}.& \left(\mathrm{S12}\right)\end{array}$

K.sub.Δd in Eq. (S17) can be evaluated from the model parameters (see next section). Nevertheless, for simplicity we take it to be equal to the corresponding out-of-plane force constant in graphite K.sub.Δd˜5 N/m.sup.44.

[0165] The resulting relative displacement is plotted in FIG. 7E versus ε/q.sup.2. Crucially, it changes sign when

[00024]${\epsilon }_{\mathrm{eff}}=\frac{\epsilon }{q^2}\sim 3.5m\mathrm{eV}.$

The experimentally measured voltages indicate that Δd in AB stacked bilayer h-BN is of the order of 10.sup.−3 Å and our DFT calculations indicate that it is positive suggesting that

[00025]$\frac{\epsilon }{q^2}\sim 3\mathrm{meV}$

see FIG. 7E), similar to expected values. Notably, other layered materials, which possess different effective ε and or q values may show different quantitative and even qualitative polarization.

[0166] Finally, it should be noted that our simplistic classical approach is sufficiently flexible to allow to study additional effects such as the dependence on the number of layers as well as external perturbations like pressure or electric field, as well as an additional in-plane component of the polarization, which we leave for future work.

Analytic Estimate of the Normal Spring Constant

[0167] To determine

[00026]$K_{\Delta d}=\frac{d^{2}E}{{d\left(\Delta d\right)}^2},$

we note that it has contributions from the interlayer LJ and Coulomb forces, as well as from the intra-layer forces. Its interlayer LJ contribution is obtained from the corresponding contribution to the energy

[00027]$E^{LJ}\left(\Delta d\right)=\frac{1}{2}\left[V_{11}^{LJ}\left(h\right)+V_{12}^{LJ}\left(h-2\Delta d\right)\right]+\frac{1}{2}\left[V_{12}^{LJ}\left(h\right)+V_{12}^{LJ}\left(h+2\Delta d\right)\right].$

The first (second) term represents the interaction of atoms at eclipsed (hollow) sites, and

[00028]$V_{11\left(12\right)}^{LJ}\left(h\right)=4\epsilon {.Math.}_{{\stackrel{.fwdarw.}{R}}_{11\left(12\right)}}\left(\frac{{\sigma }^{12}}{{\left({\stackrel{.fwdarw.}{R}}_{11\left(12\right)}^2+h^2\right)}^6}-\frac{{\sigma }^6}{{\left({\stackrel{.fwdarw.}{R}}_{11\left(12\right)}^2+h^2\right)}^3}\right).$

One can obtain a good approximation for the spring constant

[00029]${\left(\frac{d^{2}E}{{d\left(\Delta d\right)}^2}\right)}_{\Delta d=0}$

by replacing the sums Σ.sub.{right arrow over (R)}.sub.11(12) by integrals, as in F.sub.continuum.sup.LJ(h). This procedure yields

[00030]$\begin{array}{cc}{K}_{\Delta d}^{\mathrm{LJ}}=\frac{64\pi }{\sqrt{3}{R}_0^2}\epsilon .& \left(\mathrm{S13}\right)\end{array}$

[0168] For ε=3 meV this yields a spring constant of K.sub.Δd.sup.LJ=3 N/m. While additional contributions are expected from intra-layer interactions, as well as Coulombic inter-layer interactions, this value is comparable to the measured out-of-plane force constant in graphite K.sub.Δd˜5 N/m.sup.44. Since Eq. S13 captures the order of magnitude of K.sub.Δd, we can use

[00031]$K_{\Delta d}=\frac{64\pi }{\sqrt{3}{R}_0^2}\epsilon .Math.a,$

with a factor of order unity

[00032]$a\sim \frac{5}{3},$

and obtain an approximate expression for the relative displacement, see Eq. (S17), fully in terms of our model's parameters,

[00033]$\Delta d\sim \frac{F_{\Delta d}}{K_{\Delta d}}=\frac{\sqrt{3}{R}_0^2}{32\pi }\frac{\left(F_{11}^C\left(h\right)-F_{12}^C\left(h\right)\right)+\left(F_{12}^{LJ}\left(h\right)-F_{11}^{LJ}\left(h\right)\right)}{4\epsilon }.$

[0169] Having thus described several embodiments for practicing the inventive method, its advantages and objectives can be easily understood. Variations from the description above may and can be made by one skilled in the art without departing from the scope of the invention.

[0170] Accordingly, this invention is not to be limited by the embodiments as described, which are given by way of example only and not by way of limitation.