Charge-polarized interfacial superlattices in marginally twisted hexagonal boron nitride

When two-dimensional crystals are brought into close proximity, their interaction results in reconstruction of electronic spectrum and crystal structure. Such reconstruction strongly depends on the twist angle between the crystals, which has received growing attention due to interesting electronic and optical properties that arise in graphene and transitional metal dichalcogenides. Here we study two insulating crystals of hexagonal boron nitride stacked at small twist angle. Using electrostatic force microscopy, we observe ferroelectric-like domains arranged in triangular superlattices with a large surface potential. The observation is attributed to interfacial elastic deformations that result in out-of-plane dipoles formed by pairs of boron and nitrogen atoms belonging to opposite interfacial surfaces. This creates a bilayer-thick ferroelectric with oppositely polarized (BN and NB) dipoles in neighbouring domains, in agreement with our modeling. These findings open up possibilities for designing van der Waals heterostructures and offer an alternative probe to study moiré-superlattice electrostatic potentials.


Supplementary Note 2. Electrostatic imaging
The triangular potential modulation was detected using electrostatic force microscopy (EFM) and Kelvin probe force microscopy (KPFM). Both techniques are non-contact scanning probe techniques that probe the electrostatic interaction between a conductive AFM tip and the sample 3,4 . EFM measures local electrostatic force variations that can originate from either a variation in the surface potential or in the dielectric properties of the sample. In this section, we present additional information and images using EFM. KPFM is just an advanced EFM mode in which surface potential variations are recorded using an additional feedback loop, and we discuss it in the next section (see S3).
EFM images in the main text were taken using dc-EFM mode (also known as phase-EFM) 5 by simply applying a dc voltage, as illustrated in Supplementary Fig. 3a. The cantilever was oscillated at its free mechanical resonance with a dc voltage applied between the tip and the silicon substrate, and the phase shift, ∆߮, of the mechanical oscillation of the cantilever was recorded while scanning the surface. The electrostatic force experienced by the cantilever can be written as ‫ܨ‬ ሺ‫ݖ‬ሻ ൌ ‫ܥ߲‬ ‫ݖ߲‬ ⁄ ሺܸ ௗ െ ܸ ௦ ሻ ଶ /2, where C is the total tip-sample capacitance, ܸ ௗ is the applied dc voltage between the AFM tip and the sample substrate, Vs is the surface potential and ‫ݖ‬ is the tipsurface distance. To a first approximation, the phase shift directly depends on the force gradient as ∆߮ሺ‫ݖ‬ሻ ൌ െ ொ డி డ . Thus, in the presence of an electrostatic force, it can be written as Q is the quality factor of the cantilever and k its spring constant. Hence, the phase shift varies proportionally to the square of the tip-surface potential difference and to the second derivative of the tip-surface capacitance, ߲ ଶ ‫ܥ‬ ‫ݖ߲‬ ଶ ⁄ . The latter is a complex function of the geometric and dielectric properties of the tipsample system.

Supplementary Figure 3 | Schematics of the EFM setup. (a)
In dc-EFM, the AFM phase was recorded during the second pass with a dc voltage applied between the AFM probe and the silicon substrate. (b) In ac-EFM, an ac voltage of frequency ߱ was applied between the AFM probe and the silicon substrate during the second pass. The amplitude of ∆߮ሺ߱), which depends on the surface potential, showed periodic triangular modulation on the twisted-hBN crystals, while the amplitude of ∆߮ሺ2߱ሻ, which depends only on the surface dielectric properties, showed no periodic pattern.
We recorded ∆ϕ using the standard two-pass method. First, we acquired the topography image with no applied voltage. Then we retraced it with an applied dc bias of 2-3 V, the AFM feedback control switched off and the tip lifted up a few nm with respect to the first pass. The lift height, zlift, set in the range 3-5 nm, was chosen large enough to avoid short-range forces, but as small as possible to minimize the tip-surface distance, ‫.ݖ‬ This is important to maximize the lateral resolution of the technique, which is set by the ߲ ଶ ‫ܥ‬ ‫ݖ߲‬ ଶ ⁄ term and decreases with the tip moving away from the surface. By keeping the oscillation amplitude in the range of 5-10 nm, we typically took the EFM images at a total scan height of 8-20 nm from the interface between the twisted hBN crystals. We note that by probing the force gradient through ∆߮ instead of the force as in standard amplitudemodulation EFM or KPFM, the dc-EFM mode employed here is less sensitive to long-range forces from the tip cone and cantilever. Thus, it allows higher lateral resolution and it is advantageous here to study small domains. To discriminate whether the observed triangular pattern was a built-in potential or a change in the dielectric properties of the heterostructure, we took ∆ϕ curves as a function of the applied dc bias in the centre of the triangular domains, as detailed in section S3. They show the expected parabolic behaviour (see Supplementary  Fig. 7) [6][7][8] . While the curvature, set by the capacitive term ߲ ଶ ‫ܥ‬ ‫ݖ߲‬ ଶ ⁄ , was independent of the domain, the maximum of the parabola shifted a few hundreds of mV with the domain polarity. This allowed us to conclude that the triangular modulation detected in the EFM images is a built-in potential due to the interfacial charge distribution between twisted hBN crystals.
To support this conclusion, we also took images in ac-EFM mode, illustrated in Supplementary Fig. 3b, which allowed us to record dielectric images. In this mode, an ac voltage bias of frequency ߱ and amplitude ܸ was applied between the AFM tip and the sample substrate. This modulates the electrostatic force at ߱ and 2߱ with amplitude ‫ܨ‬ ሺ߱ሻ = ‫ܥ߲‬ ‫ݖ߲‬ ⁄ ܸ ௦ • ܸ and ‫ܨ‬ ሺ2߱ሻ = ‫ܥ߲‬ ‫ݖ߲‬ ⁄ • ܸ ଶ /4, respectively. While the ߱ harmonic is again proportional to both the surface potential Vs and the dielectric properties through the capacitive term ‫ܥ߲‬ ‫ݖ߲‬ ⁄ , the 2ω harmonic depends only on the ‫ܥ߲‬ ‫ݖ߲‬ ⁄ term 9 . Using two additional lock-in amplifiers, we measured the amplitude of both harmonics. Note that, also in this mode, we measured the phase shift ∆߮ of the mechanical oscillation of the cantilever instead of the force, thus recording the two phase harmonics, ∆߮ሺ߱ሻ and ∆߮ሺ2߱ሻ. Again, this is advantageous because ∆߮ is proportional to ߲ ଶ ‫ܥ‬ ‫ݖ߲‬ ଶ ⁄ and therefore it allows higher spatial resolution. Supplementary Figure 4b shows a representative ac-EFM image at ߱ on one of our twisted hBN samples. The potential modulation extends over large regions with regular and irregular triangular domains, similarly as the one observed in dc-EFM in Fig. 1d-f, and only on one side of a monolayer step. Supplementary Figure 5 shows zoom-in ߱ and 2߱ images in flat regions around bubbles filled with contamination. The domains are visible in the ߱ image, but not in the 2߱ image which depends only on the surface dielectric properties. This confirms that the triangular domains reflect a built-in potential that originates at the interface between twisted hBN crystals, not a change in dielectric properties. Figure S6 shows additional images taken in regions around other monolayer steps. Again, they show small triangular domains only on one side of monolayer steps.

Supplementary Note 3. Experimental quantification of potential modulation and KPFM imaging
We quantified the triangular potential modulation, ∆Vs, in Fig. 3f by measuring the AFM phase shift, ∆߮, as a function of the applied dc bias in the centre of two neighbouring domains, as shown in Supplementary Fig. 7. The observed ∆߮ሺܸሻ curves ( Supplementary Fig. 7c), which were acquired at the same scan height, are parabolic with same curvature set by the capacitive term ߲ ଶ ‫ܥ‬ ‫ݖ߲‬ ଶ ⁄ . On the other hand, the maximum of the parabola shifted horizontally with the domain polarity, indicating a change in surface potential [6][7][8] . We thus quantified ∆Vs as the difference between the maximum of the two parabolas. The horizontal shift of the parabola also explains the contrast inversion upon changing the sign of the dc bias ( Supplementary Fig. 7a and b). We found ∆Vs = 240 ± 30 mV on all our samples. The value was robust against variations in the domain shape, orientation and size, from micrometre range down to ∼ 30 nm, the smallest domain we could detect within our resolution. For largest domains (1 µm range) with irregular shape, the domain size in Fig. 3f indicates the smallest side of the domain.
The extracted ∆Vs was also independent of the number of layers in the hBN crystals (within the range of thickness studied here). Importantly, we confirmed the observed value of ∆Vs across different samples and using different AFM probes, thus proving the generality of our observation. We note that the maximum of the parabolas can shift with changing the sample, which in turn makes the quantification of the surface potential, Vs, difficult. However, what we quantified here and compared between samples is the variation of the potential, ∆Vs, not its absolute value, by measuring the shift between pairs of parabolas taken on the same sample. Such variation is robust against experimental conditions of different samples. In particular, it is independent of the cleanness of the surface, which directly affects Vs. We found appreciable changes in ∆Vs only when we brought the tip into contact or near contact with the surface, causing charge injection. Thus, the data in Fig. 3f were taken by carefully avoiding tip-surface contact.
It is important to note that the EFM contrast is not expected to depend on the moiré size using the experimental conditions used here, and this is confirmed by the experimental data in Fig. 3f. This is justified by the fact that Supplementary Figure 7 | Experimental potential variation using dc-EFM. (a,b) Representative dc-EFM images on marginally twisted hBN with -1,5 V and +1.5 V dc bias applied between the tip and the sample substrate. The triangular contrast reverses upon changing the sign of the dc bias. (c) dc-EFM signal, ∆߮, as a function of the applied dc bias in the centre of two neighbouring domains in (a) and (b), marked as 1 and 2. The two curves show a parabolic behaviour with same curvature, but they are shifted in the horizontal direction, which indicates a variation in the surface potential. The dashed lines are a guide to the eye.
we took the EFM images by scanning a sharp AFM tip (radius R ≈ 5 nm, cone angle θ ≈ 10º) very near to the domains plane (scan height hs ≈ 10 nm). In these experimental conditions, the EFM contrast originates in the short-range interaction of the tip with the surface area just below and around the apex. At the same time, the domains that we analysed here were much larger than or approximately equal to 30 nm in lateral size, that is, much larger than the characteristic size of the tip. Therefore, the field lines were not sensitive to the lateral extension and shape of domains larger than few tens of nm (see Ref. 10 ). It is important to note that the EFM contrast can be size-dependent on the smallest domains (30 nm < L < 100 nm) depending on the scan height, hs, set during the EFM image. Therefore, the value of hs needs to be chosen carefully. This is because on the smallest domains the built-in potential decays rapidly for values of hs typically used in EFM/KPFM (tens of nm). We have experimentally determined how the EFM contrast decays when the distance between the tip and the sample become comparable with the size of the domains. Supplementary Figure 8  carefully measured all the domains at approximately hs ≈ 10 nm height (light-red coloured region in Supplementary Fig. 8a), where ∆Vs is independent on the lateral size of the domains. This has guaranteed that we could detect the same potential difference for different moiré sizes down to ∼ 30 nm.
To corroborate the observed value of ∆Vs, we also quantified its value by KPFM imaging 4 . To do that, we used the ac-EFM setup described above (Supplementary Fig. 3b) with an additional feedback loop and a dc bias between the tip and the silicon substrate. While the tip is scanning, the dc bias is continuously adjusted by the feedback to nullify the amplitude of the ω harmonic, now equal to ‫ܨ‬ ሺ߱ሻ ൌ ‫ܥ߲‬ ‫ݖ߲‬ ⁄ ሺܸ ௗ െ ܸ ௦ ሻ • ܸ . The KPFM image thus yields the surface potential Vs of the sample, mapping its variation across all the domains, not only in their centre. We note that also in KPFM we recorded the phase shift ∆߮ሺ߱ሻ instead of the force to increase our lateral resolution. Supplementary Figure S9 shows representative KPFM images of the twisted-hBN sample shown in Fig. 1c, taken at the same scan height as in Fig. 3f (∼ 9 nm). We found the same periodic pattern in KFPM images as in the EFM images, with large areas of regular ( Supplementary Fig. 9b) and irregular ( Supplementary Fig. 9a) in (a) and (b), respectively. The surface potential variation between the centres of two neighbouring domains agrees with the value obtained in Fig. 3f from ∆߮ሺܸሻ curves. Acquisition parameters: oscillation amplitude 6 nm; lift height zlift = 3 nm; ac voltage bias of 4 V at 7 kHz.

Supplementary Note 4. Theoretical calculations
Theoretical results, with an additional analysis of the band structure, will be presented in Ref. 11 Here we shall just present the relevant results from that work.
Relaxation. The theoretical calculations were performed in two stages: first, we use LAMMPS 12 to minimize the energy using a classical potential model for relaxation 13 , using the 'inter-layer potential'(ILP) from Refs. 14,15 with the Tersoff in-layer potential 16,17 . We minimise the positions for a supercell commensurate with the hBN one, keeping the size of the supercell fixed. Alignments are plotted using an extension of the method in Ref. 13 where we take into account all six alignment options (see Fig. 2 in the main text). This leads to a strain in all these systems that is concentrated along the zone boundaries, and gives rise to a piezoelectric charge, as shown in Supplementary Fig. 10. Tight-binding model. Using the deformed positions, we then perform a tight-binding model. We neglect the modification of the in-layer hoppings due to the small bond stretching, and use a constant in-layer nearest neighbour hopping ‫ݐ‬ ൌ 2.33 eV 18 . We use a simple electronic coupling using an exponential Koster-Slater interlayer model, where ܺ and ܻ label the atomic species, ݀ ൌ 0.333 nm is the interlayer distance, and the inverse range ߙ ൌ 44 nm ିଵ .
We then diagonalize the resulting tight-binding model, either for energies near the gap (which allows us to use sparse matrix techniques, and thus study much larger moirés) or by finding all states, which is required to describe the charge density. This is calculated by summing over all occupied states, which limits the smallest angle we can perform calculations for to about 1°. Further details of these calculations, as well as further results on the electronic structure of hBN can be found in Ref. 11 .
We note that in contrast to the experimental results, all our theoretical results are for a bilayer system. However, we have shown 19 that relaxation in multilayer systems still shows a sizeable reconstruction at the interface. Nevertheless, we expect that the results we get for the piezoelectric charge may be a substantial overestimate of their real magnitude.