Two-dimensional ferroelectricity in a single-element bismuth monolayer

Ferroelectric materials are fascinating for their non-volatile switchable electric polarizations induced by the spontaneous inversion-symmetry breaking. However, in all of the conventional ferroelectric compounds, at least two constituent ions are required to support the polarization switching1,2. Here, we report the observation of a single-element ferroelectric state in a black phosphorus-like bismuth layer3, in which the ordered charge transfer and the regular atom distortion between sublattices happen simultaneously. Instead of a homogenous orbital configuration that ordinarily occurs in elementary substances, we found the Bi atoms in a black phosphorous-like Bi monolayer maintain a weak and anisotropic sp orbital hybridization, giving rise to the inversion-symmetry-broken buckled structure accompanied with charge redistribution in the unit cell. As a result, the in-plane electric polarization emerges in the Bi monolayer. Using the in-plane electric field produced by scanning probe microscopy, ferroelectric switching is further visualized experimentally. Owing to the conjugative locking between the charge transfer and atom displacement, we also observe the anomalous electric potential profile at the 180° tail-to-tail domain wall induced by competition between the electronic structure and electric polarization. This emergent single-element ferroelectricity broadens the mechanism of ferroelectrics and may enrich the applications of ferroelectronics in the future.

Ferroelectric materials are fascinating for their non-volatile switchable electric polarizations induced by the spontaneous inversion-symmetry breaking. However, in all of the conventional ferroelectric compounds, at least two constituent ions are required to support the polarization switching 1,2 . Here, we report the observation of a single-element ferroelectric state in a black phosphorus-like bismuth layer 3 , in which the ordered charge transfer and the regular atom distortion between sublattices happen simultaneously. Instead of a homogenous orbital configuration that ordinarily occurs in elementary substances, we found the Bi atoms in a black phosphorous-like Bi monolayer maintain a weak and anisotropic sp orbital hybridization, giving rise to the inversion-symmetry-broken buckled structure accompanied with charge redistribution in the unit cell. As a result, the in-plane electric polarization emerges in the Bi monolayer. Using the in-plane electric field produced by scanning probe microscopy, ferroelectric switching is further visualized experimentally. Owing to the conjugative locking between the charge transfer and atom displacement, we also observe the anomalous electric potential profile at the 180° tail-to-tail domain wall induced by competition between the electronic structure and electric polarization. This emergent single-element ferroelectricity broadens the mechanism of ferroelectrics and may enrich the applications of ferroelectronics in the future.
Ferroelectrics are well known for their applications in the non-volatile memories 4 and electric sensors 5 , and their applications have been extended to the areas of ferroelectric photovoltaics for the efficient renewable energy harvesting 6 and synaptic devices for the powerful neuromorphic computing 7 . Recently, research on ferroelectrics has been expanded to the two-dimensional (2D) limits with distinct performance [8][9][10] , including perovskite ferroelectrics at the unit-cell thickness 11,12 , single layer ferroelectrics with in-plane or out-of-plane polarization 13,14 and 2D moiré ferroelectrics by the van der Waals stacking 15,16 .
Normally, ferroelectric materials are compounds that consist of two or more different constituent elements 1,2 . The electron redistribution during the chemical bond formation instantaneously renormalizes the valence orbitals and yields the anion and cation centres. Further relative distortion, sliding or charge transfer between the positive and negative charge centres in a unit cell produces the ordering of electric dipoles to sustain the ferroelectricity [17][18][19] . By contrast, as the atoms in a unit cell of an elementary substance are identical, ordered electric dipole or even ferroelectric polarization seem difficult to form spontaneously. The realization of single-element ferroelectricity also lacks experimental demonstration so far. Nevertheless, elements situated between metals and insulators in the periodic table show flexible bonding abilities to adopt several states in one system, such as Sn atoms in the 2D Sn 2 Bi honeycomb structure show binary states 20,21 . Even in elemental boron, ionicity with inter-sublattice charge transfer was found to arise from the different bonding configurations in each sublattice (B 12 and B 2 ) 22 . The subtle balance between metallic and insulating states in these elements is easy to shift by the different sublattice environments so that both states may be realized simultaneously in a unit cell, providing possibilities to produce cations and anions in a unit cell to achieve ferroelectricity in single-element materials. Recently, some theoretical works have been devoted to exploring single-element polarity or ferroelectricity in elemental Si (ref. 23), P (ref. 24,25), As (ref. 25), Sb (ref. 25,26), Te (ref. 27) and Bi (ref. 25,26). In particular, Xiao et al. predicted that the family of group-V single-element materials in a 2D van der Waals form 28 , that is, the monolayer As, Sb and Bi in the anisotropic α-phase structure, have a non-centrosymmetric ground state to support both cations and anions in a unit cell and produce in-plane ferroelectric polarization along the armchair direction.
Article weak hybridization between the 6s and 6p orbitals so that it features partial sp 2 character other than the homogenous tetrahedral sp 3 configuration that exists in the black phosphorus 26,28 . This brings in a small buckling (Δh) between the neighbouring sublattices with the loss of centrosymmetry 30,31 . As shown in Fig. 1a, the breaking of inversion symmetry allows BP-Bi to adopt two domain states, either the Δh = d 0 or Δh = −d 0 state. Our first-principles calculations reveal the two states can be switched to each other by crossing a small energy barrier of 43 meV per unit cell (Fig. 1b). Moreover, the buckling degree of freedom enlarges the band gap and lifts the degeneracy of the p z orbitals at sublattice A and B (Extended Data Fig. 1). Taking Δh = d 0 for instance (Fig. 1c, Δh adopts the right minimum of the double-well potential), the valence band and conduction band at the Γ point are mainly contributed by the p z orbitals of A and B sublattice, respectively. When the Fermi level crosses the band gap, the valence p z orbitals at the A sublattice are fully occupied, and the p z orbitals at the B lattice are empty. In real space, this corresponds to electron transfer from sublattice B to sublattice A, leading to a spontaneously polarized character (Fig. 1d). The anharmonic double-well potential with a small barrier implies the possible ferroelectric switch between the two domain states. The nearly linear dependence between the polarization (P) and the buckling degree (Δh), along with the mirror (glide) symmetry refer to the (01) surface (in-plane central surface) (Fig. 1a), indicate that Δh is an order parameter to characterize the ferroelectricity of BP-Bi.
Experimentally, we grew BP-Bi on highly oriented pyrolytic graphite (HOPG). Figure 1e shows a scanning tunnelling microscopy (STM) image of the monolayer BP-Bi with some second layer nanoribbons along the [01] direction (Extended Data Fig. 2a) 32 . The atomic-resolved noncontact atomic force microscopy (AFM) measurement indicates two different states in two neighbouring domains separated by a domain wall (Fig. 1f,g and Extended Data Figs. 3f-o and 8). We performed force spectroscopy measurements (the Δf(z) spectra) to find out the relative distance between the Bi atoms and tip apex 33 . At the constant-height mode, the measured height difference of the turning points of the Δf(z) spectra on sublattice A and B in Fig. 1i quantitatively present the buckling degree, and the result Δh 0 = −Δh 1 = 40 pm is determined (Extended Data Fig. 5a-d). The dI/dV spectrum at the domain wall shows the p z bands (E i and E ii ) moving to a higher binding energy compared to the normal domain position, and a sharp peak occurs in the band gap (Fig. 1h), which will be elaborated in detail later.

In-plane polarization
In the buckling structure, the charge transfer between the sublattice A and B of BP-Bi is represented by the energy splitting of the p z orbitals and the variation of the surface electrostatic potential between A and B. According to the calculated band structure in Fig. 1c, the dI/dV peak corresponding to p z orbitals should be below zero (occupied state, E i ) for one sublattice but above zero (empty state, E ii ) for the other sublattice. Figure 2b shows the dI/dV cascade along the red dashed arrow (across an ABABA lattice) in Fig. 2a, revealing two traces of peaks. The valence band peak E i is clearly strong at the A sublattice (that is, points 1, 3, 5) but weak at the B sublattice (points 2, 4), whereas the conduction band peak E ii shows the opposite behaviour. At the 2D scale, the dI/dV mapping of the valence band (Fig. 2c) and conduction band (Fig. 2d) show the same feature: filled p z orbitals localize at the A sublattice, while the empty p z orbitals localize at the B sublattice, confirming the predicted electron transfer from B to A.
- [10] [01] - Another piece of evidence to prove electron transfer is the disparity in the local contact potential difference (LCPD) on each sublattice at atomic scale 34,35 . Figure 2e presents a typical Kelvin probe force microscopy (KPFM) measurement implemented by recording the frequency shift as a function of the sample bias voltage. The electrostatic force caused by the contact potential difference (CPD) between the sample and tip changes by sweeping the bias voltage. When the CPD is totally compensated by the bias (V CPD = V * ), the electrostatic force reaches the minimum, corresponding to the parabolic apex in the frequency shift curve in Fig. 2e. By recording the bias voltage V * site by site in a lateral grid, the surface potential is mapped out and shown in Fig. 2g. AFM was performed simultaneously to determine the atomic structure at the same area (Fig. 2f). It is obvious that the two sublattices A and B have different surface potentials. Using the spontaneous buckled model (Fig. 1a), we calculated the electrostatic potential shown in Fig. 2h, which reproduces the experimental LCPD map and indicates the electron enrichment at the topmost A sublattice. Ultimately, on the basis of the above dI/dV and LCPD measurement combined with the in-plane distorted atom structure, the in-plane polarization can be confirmed.

Ferroelectric switching
The polarization of a ferroelectric material can be reversed by an applied electric field 36 . In our experiments, we use the in-plane component of electric field from the conductive STM/AFM tip to switch the polarization of small ferroelectric domains close to the tip 37 . To facilitate the switching, we set the tip near a domain wall so that the size of the switched domain is comparable to the limited range of the electric field produced by the tip end (Fig. 3a). During the sample bias sweeping at a specific tip height, the IV spectra are recorded as shown in Fig. 3d. When the voltage ramps from negative to positive bias, a low conductance appears at around 0.2 V due to the band gap of BP-Bi under the tip (blue series curves). Continuous bias increase causes a small current jump labelled by the red vertical bars. While sweeping the voltage back (from positive to negative), a larger hysteretic current is maintained with a substantial conductance emerging at around 0.2 V in the gap (red series curves). According to the dI/dV curves of the domain wall ( Fig. 1h), the large in-gap current indicates the domain wall has been moved to the tip position during the application of the positive bias. When the bias voltage reaches a negative value V SW , the current jumps to the original level, suggesting the domain wall is moved back. The AFM images (Fig. 3b,c) after the forwards and backwards voltage sweeping demonstrate the movement of the domain wall directly. Figure 3c shows the domain wall movement to the left side with a distance of four lattice periods after the forwards manipulation, and Fig. 3b shows the domain wall moving back to the original position after changing the bias backwards to below V SW . Accordingly, the ferroelectric switching between the two domain states (Δh = d 0 and Δh = −d 0 ) is experimentally verified. The domain wall in Fig. 3a can be assigned to the 180° head-to-head ferroelectric domain wall with the polarization labelled in Fig. 3b,c.
In the hysteretic domain manipulation, the switching voltage at both positive bias side (V SW0 ) and negative bias side (V SW ) show tip-height dependence. When the relative tip height Δz rises from 0 to 60 pm, the switching voltage V SW0 at the positive bias side is found to increase while V SW at the negative bias decreases (Fig. 3d,e). We also measured the LCPD at different heights (V CPD ) (Fig. 3e), which reveals a gradual LCPD decrease with the tip height (Δz) climb due to the semimetallicity of graphite 38 . The tip in our setup is grounded with a bias voltage (V S ) applied to the sample, the electric potential of the HOPG can be written as Φ S = V S − V CPD . Thus, the increase of tip height strengthens the electric field under the tip through the decrease of the V CPD , but weakens the electric field by the rise of the tip-sample distance. Nonetheless,  Article because Φ S is small at the negative bias side and large at the positive bias side, we find the electric field is dominated by the LCPD at the negative bias side (V SW ), and dominated by the tip-sample distance at the positive bias side (Extended Data Fig. 7g,h). Therefore, the switching bias required to trigger the same domain switch at a higher tip height has to decrease at the negative bias side but increasing at the positive bias side accordingly, which also explains the correlated variation between V SW and V CPD in Fig. 3e.
Quantitatively, the electric potential and electric field below the tip are derived by solving the Poisson equation in prolate spheroidal coordinates. The switch bias at the positive bias side was set to V S = 1.0 V (approximately Φ S = 1.6 eV). The calculated results are shown in Fig. 3g. Differentiating the electric potential on the BP-Bi surface (Φ) produces an in-plane electric field of as large as −20 mV Å −1 at the distance of roughly 10 Å to the centre (r = 0 Å). On the other hand, we can estimate the coercive field by E α β = 2 − /27 The estimated E c = 15.7 mV Å −1 is comparable to the calculated in-plane electric field above (Fig. 3g). This may be understood that the highly localized in-plane electric field works mainly on the domain part between tip and wall, thus demanding that the electric field overcomes E c . It is noticeable that they are also close to that in the reported ferroelectric polymer PVDF 39 and strained transition metal oxides 40,41 . As spontaneous polarization (P s ) is determined by the same α and β through P α β = − / s , the consistency between E c and in-plane electric field E in our experiment also verifies the DFT result P s = 0.41 × 10 −10 C per m (Fig. 1b), which resembles the polarization of monolayer SnTe (roughly 1 × 10 −10 C per m) 42 .

180° domain walls
Besides the 180° head-to-head domain wall in the monolayer BP-Bi (Figs. 1e, 3a, 4a and Extended Data Fig. 9), we also observe the conjugated 180° tail-to-tail domain wall (Fig. 4c). The AFM measurements show that the 180° head-to-head domain wall has a wall width of about 15 Å (three unit cells, shadow area in Fig. 4e), whereas the 180° tail-to-tail domain wall has a wider width of roughly 56 Å (12 unit cells, shadow area in Fig. 4g). The dI/dV measurements of the head-to-head domain wall reveal a one-dimensional electronic state in the band gap, and both conduction band and valence band near the wall bend to a higher binding energy (Fig. 4b), indicating the accumulation of electrons around the wall. A close determination of the band bending by extracting E i peak and measuring the LCPD by KPFM show the same bending profile and a total band bending of about 170 meV (Fig. 4f). When turning to the tail-to-tail domain wall, the dI/dV spectra, however, show the incongruous band movement in the valence band and conduction band (Fig. 4d). The same band bending assessment methods by the E i analysis and LCPD measurement in Fig. 4h suggest a smaller but identical bending direction to the head-to-head domain wall (Fig. 4f).
To understand the experimental observation above, we calculate the order parameter profile at the domain wall using the Landau-Ginzburg-Devonshire theory 43,44 . Considering the intrinsic p-type doping (Extended Data Fig. 4), the screened Coulomb interaction is included to take into account the screening effect of charge carriers on the ferroelectric instability 45 . Figure 4i-l shows the calculated wall width and surface potential profile in the 180° head-to-head and tail-to-tail domain walls. It is apparent that there is nearly no difference between the cases with or without the screened Coulomb interaction in the head-to-head domain wall (Fig. 4i,j). The wall width is about three unit cells and the surface potential increase exponentially at the domain wall, consistent with the experimental observations (Fig. 4e,f). However, for the tail-to-tail domain wall, the wall width and the surface potential have totally different behaviours regarding the cases with and without a screened Coulomb interaction (Fig. 4k,l). For the case without a screened Coulomb interaction, the tail-to-tail domain wall has a similar narrow width (three unit cells) as the head-to-head domain wall. Involving the screened Coulomb interaction results in a wall width four times larger (12 unit cells), which matches the experimental observations.
As demonstrated, BP-Bi is heavily p-type doped (Extended Data Fig. 4). The upwards band bending at the tail-to-tail domain wall in principle will increase the local carrier concentration drastically, which accordingly decrease the Thomas-Fermi screening length and strongly screen the Coulomb interaction or depress the spontaneous polarization P s (refs. 46,47). Generally, the wall width can be simplified as D k P β = 4 /2 W s 2 (k is the gradient coefficient) 36 , thus the rise of charge carrier density at the domain wall broadens the wall width and flattens the surface potential at tail-to-tail domain wall (red dashed lines in Fig. 4k,l). By contrast, the band bending at the head-to-head domain wall makes the Fermi level shift to the middle of the band gap and lessens the carrier density, which consequently narrows the wall width (Fig. 4i). However, as the screening length at low carrier density is notably larger than the efficient range of the Coulomb interaction 45 , the screening effect changes little with the continued reduction of the carrier concentration, resulting in a slight width shrinking at the headto-head domain wall.
In addition, the evolution of the electronic structure correlated to the atomic buckling (Δh) at the tail-to-tail domain wall is discussed. Early theoretical and experimental research showed that the band structure is highly buckling dependent 30,48 . Our refined DFT calculations suggest that the position of the Fermi level decreases monotonously with the buckling reduction (Extended Data Fig. 1c), which implies an increase in the work function or gradual charge transfer between the substrate and BP-Bi layer at the wide tail-to-tail domain wall. With this consideration, the amended calculation of the surface potential reverses at the tail-to-tail domain wall, and acts as a similar band bending as that of the head-to-head domain wall (red solid curve in Fig. 4l). This result reproduces the electric potential measurements (Fig. 4h) and Article illustrates the conjugative correlations between electronic structure and ferroelectric distortion in this monolayer BP-Bi system.

Conclusion
In this work, we experimentally confirmed in-plane polarization and observed ferroelectric switching in a single-element BP-Bi monolayer, demonstrating the capability to realize ferroelectric polarization in an elementary substance or single-element compound 22 (for example, bismuth bismuthide here). The spontaneous charge redistribution and the ferroelectric atomic distortion in the BP-Bi show the intercorrelation between the electronic structure and inversion symmetry.
Owing to the heavily p-type doping of BP-Bi, the screened Coulomb interaction from the carrier density and the electronic modulation by the ferroelectric distortion were observed at the 180° tail-to-tail domain wall. The single-element ferroelectricity inspires the advantages of ferroelectrics to electrically modulate the band structure, as well as other potential emergent phenomena such as topology and superconductivity beyond magnetism 19 .

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-023-05848-5.

Sample preparation and characterization
Single layer BP-Bi was prepared in an ultra-high vacuum molecular-beam epitaxy chamber that connected to the Omicron LT-STM system (1 × 10 −10 mbar). HOPG was first cleaved in air and then immediately loaded into the molecular-beam epitaxy chamber to degas the surface contamination. Bismuth with a purity of 99.999% was evaporated by a K-cell to the substrate that kept at room temperature. All the STM or scanning tunnelling spectroscopy (STS) and AFM/KPFM measurements were performed at 4.3 K with a tungsten tip glued to the qPlus AFM sensor, except the variable temperature experiment reached the respective temperature by the liquid nitrogen cryostat cooling and resistance counter heating. The clean and atomically sharp tip was obtained by repeat voltage pulses and controllable indentations on the Au(111) substrate. Further vertical CO molecule manipulation was carried out on the same surface to get a CO functionalized apex 49 . Differential conductance (dI/dV) was obtained by a lock-in technique with a 3-10-mV and 963-Hz modulation superimposed on the sample bias. The AFM sensor with a resonance frequency of 27 kHz and a quality factor of 30,000-70,000 was excited to an amplitude of 40 pm for the AFM measurement and 60 pm for the KPFM measurement. All the images have been processed and rendered using MATLAB and WSxM software 50 .

First-principles calculations
First-principle calculations were carried out within the DFT formalism, as implemented in the Vienna Ab initio Simulation Package 51 . The local density approximation was used for exchange-correlation functional with a projector augmented-wave 52 pseudopotentials, treating Bi 5d 6s 6p as valence electrons. The cut-off of the plane wave basis was set above 500 eV to guarantee that the absolute energies are converged to 1 meV. Grimme's method 53 was used to incorporate the effects of van der Waals interactions in all geometry optimizations. The vacuum region was set to larger than 25 Å to minimize artificial interactions between images. The 2D Brillouin zone was sampled by a Monkhorst-Pack 54 k-point mesh and it is 12 × 12 for the unit cell. The positions of atoms were fully relaxed until the Hellmann-Feynman force on each atom was less than 0.01 eV Å −1 . The convergence criteria for energy were set to 10 −5 eV. The electronic-structure calculations were performed within the scalar relativistic approximation, including a spin-orbit coupling effect. Hybrid functional (HSE06) 54 was also used to confirm the band gap. All energy levels were calibrated by the corresponding core levels.

Moiré pattern and domain wall
Owing to the twofold symmetry of the BP-Bi, the moiré patterns produced by stacking single layer BP-Bi on HOPG can be classified into two categories: in type 1, the zigzag BP-Bi ([01] direction) aligns at a small twist angle θ with respect to the zigzag direction of graphite; in type 2, the zigzag BP-Bi aligns at a small twist angle θ with respect to the armchair direction of graphite. Among them, the ±1.7° type 2 (θ = ±1.7°) moiré pattern is unique for the homogeneous stacking sequence between BP-Bi and graphite along the short axis of the moiré superlattice (Extended Data Fig. 2e), which therefore shows the stripe contrast in the STM measurement (Extended Data Fig. 2d). In this kind of moiré pattern, the domain wall is constrained to along the stripe direction owing to the corresponding stripe-like strain and potential distribution in the moiré pattern (Extended Data Fig. 8a and the bottom-right part of Extended Data Fig. 2a). Conversely, the type 1 moiré pattern has the gradually changing stacking sequence along both periodic directions (Extended Data Fig. 2b,c), thus allows the 180° domain wall with an incline angle to exist (Extended Data Fig. 2a,f,g). For all the domain walls, there are two critical angles that are used to describe them adequately: the angle between polarization vectors of two neighbouring domains and the smallest angle between the domain wall and the polarization vector nearby. We name the last one as the incline angle and only use this part ahead of the domain's name when the front one is 180°. But we neglect the prefix when the incline angle is 90°, because this kind of domain wall (90° inclined 180° head-to-head or tail-to-tail domain wall) is the most frequently found (Figs. 3 and 4 and Extended Data Figs. 2a and 9). Extended Data Fig. 8 shows some uncommon domain walls in the measurement, they are charged 54° inclined 180° head-to-head domain wall (Extended Data Fig. 8a), neutral 90° head-to-tail domain wall (Extended Data Fig. 8b) and charged 90° head-to-head domain wall (Extended Data Fig. 8c).

Band edge and band gap
From the dI/dV measurement on BP-Bi domain area, it is easy to identify a band gap with the valence band maximum (VBM) near the Fermi level and the conduction band minimum (CBM) at around 0.27 eV. To determine the band structure and band edge precisely, we measured the quasiparticle interference (QPI) of the VBM and CBM at the spatial domain and energy domain (Extended Data Fig. 4). The Fourier transform of the 2D STS maps at both valence band and conduction band show the intravalley scattering characterized by the q c and q v in the centre of the reciprocal space. Fourier transform-STS along the high symmetry direction (Γ-X 2 ) show the parabolic dispersion at both valence band and valence band near the Fermi surface. The parabolic fittings of the QPI yields the VBM and CBM to be at 26 ± 5 and 289 ± 5 meV, respectively, which confirms a band gap of roughly 260 meV and the intrinsic hole doping in the BP-Bi.

Buckling determination and edge reconstruction
Extended Data Fig. 5 shows the AFM measurement of a crater and two different zigzag edges in the BP-Bi. Although the crater is so small that the edges are very short, the reconstruction can be easily confirmed on the atomic scale. Further inspection of the long straight zigzag edges in a BP-Bi ribbon island also shows the same results. The reconstructed gaps and 4a superperiod are responsible for the missing of regular band bending at both zigzag edges (Supplementary Information) 55 .
The buckling (Δh) of the BP-Bi can be roughly determined by obtaining the height difference of the turning points of the Δf(z) spectra that were measured above the A and B sublattice. A more precise measurement requires the subtraction of the background force. Here we use the force spectroscopy that was obtained at the crater as the van der Waals force background (curve C). Removing this term gives a smaller buckling value (Δh = 0.38 Å) but still a different frequency shift at the turning point of the two Δf(z) spectra, which means the background force at A and B is different. This can be understood to provide the information that the nearest and next-nearest neighbours of A and B have different z coordinates and, the A and B atoms have different electrostatic charges. At the tip-sample distance in our experiment, assuming that the electrostatic charge difference is neglectable and the background force for A and B has the same form F b (z), thereof the (subtracted) net force for atom A and B can be written as F A0 (z) = F A (z) − F b (z − Δh) and F B0 (z) = F B (z) − F b (z), respectively. Because A and B are the same Bi element, curves F A0 (z) and F B0 (z) should be identical by a translation of Δh along the z direction, which can be expressed as F A0 (z − Δh) = F B0 (z). With this consideration, we derived the net buckling Δh = 0.35 Å as shown in Extended Data Fig. 5d. Nevertheless, as the derived net buckling (Δh = 0.35 Å) is close to the raw buckling (Δh = 0.40 Å) that obtained without background-force subtraction, and the difference between them fade away when reducing the buckling to zero at the domain wall, we still use the raw buckling Δh by simply calculating the height difference of the turning points of the Δf(z) spectra to depict the wall thickness in Fig. 4.

Domain manipulation
The polarization switching at another 180° head-to-head domain wall is shown in the Extended Data Fig. 7. Similar to the one in Fig. 3, the ferroelectric hysteresis during the electric field manipulation is distinguishable by the current variation at various tip-sample distances. But we note the switching voltage (V SW ) changes from wall to wall because of the localized strain and defect, which is also the reason that causes the hysteresis loop a lateral shift so that deviates from the symmetric schematic in the inset of Fig. 3d (ref. 56).
To evaluate the electric field below the tip, we did the calculation in the prolate spheroidal coordinates (η, ξ). By treating the tip surface as a hyperboloid η t and the graphite surface as an infinite metal plane η = 0, the potential in the tunnelling gap reads: Transformation to the Cartesian coordinates (r, z) gives rise to the potential distribution Φ(r, z) and the derivative in-plane electric field E(r) at the BP-Bi plane (Fig. 3f,g).
In the meantime, the tip-height-dependent electric field at both switch sides can be calculated. As shown in Extended Data Fig. 7g,h, the electric field intensity (absolute value) below the tip decreases with the tip lifting at the positive bias side (V S = 0.7 V), but increases at the negative bias side (V S = −0.4 V). This suggests the electric field intensity is mainly dominated by the tip-sample distance elongation at the positive side, while by the decline of V CPD at the negative side, which explains the upwards shift of V SW0 , V SW2 and downwards shift of V SW , V SW1 when increasing the tip-sample distance (Fig. 3d,e and Extended Data Fig. 7f).

Ferroelectric phase transition
At higher temperatures, it is challenging to perform atomic force spectroscopy to extract the ferroelectric distortion directly. But the small atomic distortion can be magnified by the moiré pattern so that it is distinguishable by the STM topographic measurement directly. Particularly, the drastic buckling reversal at the 180° head-to-head domain wall contribute to a lateral shift of the moiré superlattice. This yields a kink of the moiré lattice crossing the wall. Extended Data Fig. 9 shows the domain wall we measured at a series of different temperatures in the same area. At low temperatures, the kink produced by the 180° head-to-head domain wall is easy to find (blue line), but it smears at 165 K and disappears at 210 K. Apart from the STM topography measurement, we did the dI/dV measurement at the domain wall, and found the peak of the in-gap state also shows similar behaviour, that is, the peak weakened with the rise of the temperature and disappeared at 210 K. The gradual changes of the kink and spectra not only exclude the tip-induced domain manipulation that could move the domain wall away by the electric field, but also derive a Curie temperature of about 210 K with a possible second-order phase transition.

Calculation of the order parameter profile at the 180° head-tohead and tail-to-tail domain wall
We use a three-dimensional (3D) counterpart that contains multilayers of 180° head-to-head or tail-to-tail domain walls along the z direction to inspect the development of the domain profile. The equivalent 3D model with an adjusted interlayer distance has the merits of involving the screening effect of the substrate and simplifying the calculation to be one-dimensional. To satisfy the free energy minimization of the ferroelectrics with second-order phase transition, we have 36,43,44 αP βP k P where x is the axis perpendicular to the domain wall plane; α, β, k have the same definitions as those in the main text and β > 0, α = α 0 (T − T C ) with the constant α 0 > 0. In the meantime, the electric potential φ and polarization P across the domain wall (along the x axis) fulfil the Poisson equation where the ε is the permittivity, ρ is the charge density that mainly contributed by the band edges of BP-Bi. According to the density of states (DOS) features reflected by the dI/dV curve and the QPI measurement in the Extended Data Fig. 4, we approach the carrier concentration by using the parabolic band edge at the CBM and a non-parabolic mode at the VBM with the fitted parameters. At T = 4.3 K under the 2D limit, they can be approximated as Here, ħ is the reduced Planck's constant, E F is the Fermi energy level, e is the absolute value of the electron charge, m C (m V ) and E C (E V ) is the effective-mass and the energy of the conduction (valence) band edge, respectively. The solving of the differential equations above produces the conjugated results at the head-to-head and tail-to-tail wall (blue curves in Fig. 4i-l).
When considering the Coulomb screening of the ferroelectric dipole, we calculate the Curie temperature at different screening length or charge concentration to include the screening term 13,57 . Thus, we rewrite the coefficient α as α 0 (T − T C (ρ)). From the results of the red dashed curved in Fig. 4k, substantial wall broadening can be observed because of the weakening of the spontaneous polarization near the wall. Further consideration of the Fermi energy at different buckling level (Δh) is conducted by introducing an extra carrier-concentration term ρ′(Δh), which is approximated linearly by referring to the DFT calculations (Extended Data Fig. 1c).

Data availability
The main data supporting the finding of this work are available within this paper. Extra data are available upon reasonable request from the corresponding authors.

Code availability
The computer code used for numerical calculation and data processing are available from the corresponding authors upon reasonable request.