Main

Real-space imaging of electric fields at the nanoscale is an important aim across many emerging fields, as near-surface fields are tied to the electrical polarization or charge distribution of the underlying system. Sensitive imaging of nanoscale electric phenomena has been demonstrated by a number of techniques, most prominently by electrostatic force microscopy15 and piezoresponse force microscopy (PFM)4, along with the related techniques of Kelvin probe force microscopy16, low-energy electron microscopy5 and emerging scanning quantum technologies1,17. However, most of these techniques are limited to low temperatures or high-vacuum conditions, require thin-film samples and back electrodes and measure indirect quantities, such as piezoelectric coefficients or surface potentials. Nitrogen-vacancy (NV) centres in diamond9,13,18 provide a path to quantitatively image electric fields under ambient conditions, do not require back electrodes or applied voltages and measure a quantity directly proportional to the surface polarization.

In this work, we apply scanning NV microscopy to map static electric stray fields above surfaces with sub-100 nm resolution. Using mechanical oscillation of the tip12 to overcome the challenges of static screening19 and low coupling to electric fields10,11, we reach an excellent sensitivity of 0.24 kV cm−1 Hz−1/2, on a par with the sensitivities demonstrated for a.c. field detection in bulk diamond11. We illustrate the impact of our approach by imaging patterned domains in application-relevant6 ferroelectric thin films and by mapping the natural domain configuration in a prototypical improper ferroelectric.

Electric field sensing with NV centres relies on a local electric Stark effect. The Stark effect causes a shift in the NV spin energy levels that is measured using optically detected magnetic resonance (ODMR)10,11. The Stark effect is anisotropic and largest in the transverse plane of the NV centre, which sits perpendicular to the anisotropy axis (zNV axis). This leads to an in-plane electric field coupling that is approximately 50× larger than the out-of-plane coupling11,13. To maximize the in-plane electric field response and simultaneously suppress the response to magnetic fields, a small magnetic bias field is applied transverse to the NV anisotropy axis. Electric field detection in this configuration has previously been demonstrated in bulk diamond, where electric fields are created with external electrodes11,20,21, charged scanning probes22,23, surface band bending24 and intrinsic dopant charges25.

In our experiment, the NV centre is embedded in the tip of a diamond scanning probe (Fig. 1a). The scanning probe arrangement allows us to extend electric field sensing to image general materials systems, including ferroelectrics. We mount the diamond probe on a quartz tuning fork oscillator providing force feedback for safe approach and scanning. Owing to the tip fabrication procedure, the NV anisotropy axis is ~35° away from the scan plane (Fig. 1b).

Fig. 1: Scanning NV electrometer.
figure 1

a, The NV centre is located at the apex of the diamond tip and is oscillated in shear mode through the electrical drive of a tuning fork (not shown) while it is scanned over the surface. Laser and microwave pulses synchronized to the drive are used for signal readout. A three-axis piezo stage underneath the sample is used for positioning. b, Geometry of tip and sample with the vector orientations of the NV spin, magnetic bias field and sample electric field. Sample surface charges and screening charges on the tip are also shown. c, The in-plane electric field detection axis (purple) is determined by the in-plane direction of the magnetic bias field (grey). (x, y, z) and (xNV, yNV, zNV) denote the laboratory and NV frames of reference (see Methods for definition). d, Spin-echo pulse scheme synchronized to the tip oscillation for the a.c. measurement of the electric field gradient.

Full size image

To enable electric field (E) detection, we accurately orient an external bias field of 5−12 mT transverse to zNV (θB = 90°) with ~0.5° of uncertainty (Methods). In this bias field configuration, the spin transition frequencies ω± are linearly sensitive to electric fields,

$${\omega }_{\pm }\approx {\omega }_{\pm }^{(0)}\mp 2\uppi {k}_{\perp }{E}_{\perp }\cos (2{\varphi }_{B}+{\varphi }_{E})$$
(1)

while correcting for magnetic fields up to second order11,13. Here, \({\omega }_{\pm }^{(0)}\) are the spin resonance frequencies in the absence of the electric field (Supplementary Section 1), k⊥ = 16.5 Hz V−1 cm is the coupling constant21, E⊥ is the magnitude of the in-plane E field vector and \({\varphi }_{B}=\arctan ({B}_{{y}_{{{{\rm{NV}}}}}}/{B}_{{x}_{{{{\rm{NV}}}}}})\) and \({\varphi }_{E}=\arctan ({E}_{{y}_{{{{\rm{NV}}}}}}/{E}_{{x}_{{{{\rm{NV}}}}}})\) are the in-plane angles of the magnetic bias field and electric field vectors, respectively (Fig. 1c). Equation (1) neglects strain interactions (Methods). The angular dependence on the bias field results in a maximal frequency shift when φE = −2φB (modulo π), defining the detection axis (Fig. 1c). To polarize, manipulate and detect the NV spin state we use a combination of laser and microwave pulses, together with a single-photon counting module (Fig. 1a,d and Methods).

An important concern with static field measurements is screening of electric fields by mobile charges on the diamond tip (Fig. 1b). This issue has hindered previous attempts at implementing a scanning NV electrometer19,26. To overcome this screening, we oscillate the diamond sensor using the mechanical resonance (f ≈ 32 kHz) of the tuning fork and detect the resulting a.c. electric field (Fig. 1d). This a.c. electric signal is proportional to the electric field gradient in the oscillation direction12. In addition to alleviating static screening, the a.c. detection also improves sensitivity by at least an order of magnitude11,12,18. For the spin-echo pulse sequence shown in Fig. 1d, the field-induced coherent phase accumulation of the NV spin is \({\phi }_{\pm }=\pm 4{k}_{\perp }{E}_{{{{\rm{ac}}}}}\cos (2{\varphi }_{B}+{\varphi }_{{E}_{{{{\rm{ac}}}}}}){\sin }^{2}(\uppi f\tau /2)/f\) where \({E}_{{{{\rm{ac}}}}}={x}_{{{{\rm{osc}}}}}\sqrt{{({\partial }_{x}{E}_{{x}_{{{{\rm{NV}}}}}})}^{2}+{({\partial }_{x}{E}_{{y}_{{{{\rm{NV}}}}}})}^{2}}\) and \({\varphi }_{{E}_{{{{\rm{ac}}}}}}=\arctan ({\partial }_{x}{E}_{{y}_{{{{\rm{NV}}}}}}/{\partial }_{x}{E}_{{x}_{{{{\rm{NV}}}}}})\), xosc is the oscillation amplitude along x and τ is the evolution time (Supplementary Section 1).

We demonstrate scanning electrometry by imaging the electric fields appearing above the surfaces of two ferroelectric materials. Our first sample is a 50-nm thick film of out-of-plane polarized lead zirconate titanate (Pb[Zr0.2Ti0.8]O3, PZT) grown on top of a SrRuO3-buffered (001)-oriented SrTiO3 substrate. PZT, the most technologically important ferroelectric, has a large polarization (P ≈ 75 μC cm−2) that is ideal for an initial demonstration. To create recognizable structures, we write a series of ferroelectric domain patterns by locally inverting the polarization of the film using a conductive atomic force microscopy (AFM) tip. Figure 2a shows an image of the out-of-plane PFM contrast corresponding to one of these patterns.

Fig. 2: Vector electrometry of a piezoelectric PZT film.
figure 2

a, Out-of-plane PFM image of a patterned eight-crossed domain structure. b, NV electrometry taken over the same region with the bias magnet oriented in-plane (B⊥zNV). c, Same as b but with the bias magnet oriented out-of-plane (B∣∣zNV). d, NV electrometry maps from a patterned square domain. The two maps show orthogonal in-plane field components, obtained by shifting the magnetic angle φB by 45°. e, Corresponding simulated electric field images with schematics of the magnetic bias angle (grey) and detection axis (purple). f, Laboratory-frame vector plot of the electric field \({{{{\bf{E}}}}}^{{{{\rm{lab}}}}}=({E}_{x}^{{{{\rm{lab}}}}},{E}_{y}^{{{{\rm{lab}}}}},{E}_{z}^{{{{\rm{lab}}}}})\) reconstructed from the maps in panel d. \({E}_{x}^{{{{\rm{lab}}}}}\) and \({E}_{y}^{{{{\rm{lab}}}}}\) components are represented as arrows and \({E}_{z}^{{{{\rm{lab}}}}}\) is shown as a colour. g, Reconstructed surface charge density σ revealing the written square domain pattern. h, Corresponding out-of-plane PFM image. Dwell times are 12 s per pixel (b,c) and 10 s per pixel (d). Scale bars, 1 μm. rad., radian; Sim., simulation; arb., arbitrary units. White dashed lines outline the inner domain.

Full size image

Figure 2b presents an NV electrometry map taken above the same location. Owing to our gradiometric detection scheme, the signal is maximum near vertical edges of the pattern. This directionality reflects the horizontal oscillation direction of the sensor, and other oscillation directions (such as a tapping mode) may be used to acquire different spatial signatures12. To verify that the signal is indeed due to electric fields, we purposely misalign the bias field to θB = 95° (Supplementary Fig. 1) and 180° (Fig. 2c). As expected, the signal disappears under the bias field misalignment as the NV centre becomes insensitive to electric fields. We additionally tried detecting the E field in a d.c. sensing mode and observed no signal (Extended Data Fig. 1). Thus, dynamic (a.c.) operation is essential for overcoming screening and enabling static electric field sensing. Detection at higher frequencies, and with multipulse measurement schemes, is shown in Extended Data Fig. 1 and Extended Data Fig. 2.

By rotating the in-plane angle φB of the bias field, we can rotate the in-plane detection axis (Fig. 1c)11. In particular, by acquiring electric field maps shifted in φB by 45°, it becomes possible to map two orthogonal components of the E field signal. Figure 2d shows such orthogonal electric gradient maps from a square domain. The experimental maps are in good agreement with numerical simulations of a square domain with constant surface charge (Fig. 2e). By combining the orthogonal field components we are able to reconstruct the full three-dimensional electric field vector above the domain (Fig. 2f and Methods). In principle, field maps such as Fig. 2f could allow measurement of the domain wall width and of its possible chiral state27. However, our spatial resolution (~100 nm; Extended Data Fig. 3) is not yet sufficient to resolve the structure of the <10 nm domain walls in PZT28. Using reverse propagation of Coulomb’s law, we also compute the equivalent surface charge density σ (from the polarization P) at the top surface of the ferroelectric film (Fig. 2g and Methods), where σ = P ⋅ n and n is the surface normal. The reconstruction involves two charge sheets of opposite sign, since further analysis (below) reveals the presence of charges beneath the surface. Our result is in excellent agreement with the out-of-plane PFM image shown in Fig. 2h, where both a defect and the asymmetric shape of the domain are reproduced.

To further interpret our measurements we develop simple models of surface charge distributions and their associated electric field gradients (Fig. 3a–d). By comparing these to the experimental line scans (Fig. 3e), we can discriminate between different surface charge scenarios and extract quantitative information on the surface charge density σ. We find that our data are best fit by the model shown in Fig. 3a, which includes two layers of opposite charges on the top and bottom of the PZT sample. The other models (Fig. 3b–d) are inconsistent with the polarity or shape of the experimental line scans. In particular, the Lorentzian shape of the \({E}_{z}^{{{{\rm{lab}}}}}\) field gradient produced by a net surface charge (models in Fig. 3b,c) fails to reproduce the dips in the signal at either side of the domain. Additionally, the models in Fig. 3c,d would result in a polarity opposite to what is expected from the known polarization of our PZT sample. The ability to distinguish between different charge models and to quantify the screening efficiency of a back electrode will be useful for analysing charge dynamics and screening behaviour at ferroelectric domain walls, interfaces and surfaces29, even once the materials are buried in a device architecture30.

Fig. 3: Surface charge models at domain walls in PZT.
figure 3

a–d, Surface charge models and negative z component of the laboratory electric field gradient (\({\partial }_{x}{E}_{z}^{{{{\rm{lab}}}}}\)) over a domain wall, which corresponds to the NV signal for φB = 121°. a, Ideal bound charge model with monopole charge sheets on the top and bottom surfaces. t is the sample thickness and for thin films the field gradient is dipole-like. b, Same as a with complete screening from the bottom electrode. c, Same as a with an adsorbed top layer separated by a distance d. The adsorbed layer screens the top surface and the signal is mainly produced by the bottom layer. d, Same as c with complete screening from the bottom electrode, resulting in a dipole surface. e, NV electrometry line scans taken across the square domain shown in Fig. 2d. Bias field angles φB are listed on the left and schematically shown (with detection axis) on the right. Profile fits (grey) are only compatible with surface charge model a. Error bars are ± shot-noise propagated uncertainties (see Methods for details). Black arrows indicate polarization direction.

Full size image

Fits to the line scans in Fig. 3e yield an effective surface charge density of approximately σ = (3.5 ± 0.3) × 10−3 μC cm−2. The fits use the known bias field direction, NV centre orientation and stand-off distance (Methods and Supplementary Section 2). This surface charge density is about four orders of magnitude lower than the expected value for our PZT sample (75 μC cm−2, ref. 31). The large difference may be attributed to a number of screening mechanisms on both the sample and diamond tip. This includes surface screening from adsorbates32 or the formation of a charged, off-stoichiometric surface layer33 at the surface of the film. On the diamond tip, screening could be a result of the dielectric constant (ϵr = 5.7) as well as partial screening from an adsorbed water layer and mobile charges still present at ~32 kHz (ref. 19). Calibration against a known electrode may allow disentangling tip and sample screening34.

We further illustrate scanning electrometry by imaging the natural domain pattern of an improper ferroelectric, the hexagonal manganite YMnO3. Hexagonal manganites are exciting benchmark materials because their surface polarization (P ~ 5.5 μC cm−2 for YMnO3) is over an order of magnitude smaller than our PZT sample and typical for ferroelectric and multifunctional materials8,14. In addition, hexagonal manganites are type I multiferroics and become antiferromagnetically ordered below TN ~ 100 K (refs. 14,35). While the antiferromagnetic domain pattern has been imaged with scanning NV magnetometry at cryogenic temperatures26, here, we focus on the ferroelectric domain pattern accessible under ambient conditions.

Figure 4a shows an electrometry map recorded above a polished, bulk YMnO3 sample. Although the signal-to-noise ratio is lower compared to the PZT data, as expected from the smaller polarization, domains (including vortex domains) are clearly visible and resemble the pattern observed by PFM (Fig. 4b). The NV electrometry image also reveals long, straight line-like features and defects, which we interpret as charge accumulation near topographic features such as polishing marks. Correlative magnetic measurements over the same region show no magnetic signal (Supplementary Fig. 2).

Fig. 4: Naturally occurring ferroelectric domain pattern in hexagonal YMnO3.
figure 4

a, NV electrometry on YMnO3. Domains appear with constant contrast, which may be attributed to tip oscillations that deviate from pure shear mode (Extended Data Fig. 1). Bright features reflect patch charges at topographic defects. b, Out-of-plane PFM on YMnO3 over the same region. Artefacts from topographic cross-talk are not observed in panel a. White arrows indicate (1) a vortex domain and (2) a 180° domain wall. Scale bar, 5 μm.

Full size image

In summary, we demonstrate nanoscale electric field imaging using a scanning NV microscope. Crucial to our experiment is the use of gradiometric detection12, which both alleviates static field screening at the diamond surface and greatly improves the sensitivity compared to static sensing schemes. The electric field sensitivity of ~0.24 kV cm−1 Hz−1/2 demonstrated in our work (Extended Data Fig. 2) is similar to that shown in the bulk11.

By imaging ferroelectric domains through their stray electric fields we demonstrate a complementary imaging method to existing scanning probe techniques. While the spatial resolution (~100 nm) demonstrated in our work is currently outmatched by these techniques4,15,16, moving the tip closer to the sample surface is expected to lower this value to below 50 nm, and perhaps below 30 nm (ref. 36). Moreover, the long measurement times of several seconds per pixel can be reduced through increased dynamical decoupling and excited-state spectroscopy at cryogenic temperatures37. The benefits of NV electrometry lie in its non-perturbing aspect, as it can quantitatively measure an unknown electric dipole configuration (electrical polarization and buried charged planes) without the need of back electrodes or applied voltages, all without topographic cross-talk in the measurement signal.

Looking towards future applications, an intriguing prospect is the correlative imaging of magnetic and electric fields in multiferroic and multifunctional materials35,38. Since the sensitivity to electric and magnetic fields can be adjusted by a simple re-orientation of the magnetic bias field, magnetic and electric field maps can be recorded from the same sample region without breaking experimental conditions. This will allow analysing multiple-order parameters in situ. Together, the multimodal capability opens exciting opportunities to study magneto-electric coupling and detail the formation and structure of domains and domain walls in these multifaceted and technologically relevant materials. Finally, the ability to probe buried charged surfaces may be instrumental in the pursuit of energy-efficient device concepts based on electric field control of magnetization39,40.

Note added in proof: While finalizing this manuscript, we became aware of related work34 describing scanning NV electrometry using a similar gradiometry technique to image electric fields from electrodes.

Methods

Experimental set-up

Experiments were performed at room temperature with a custom-built scanning NV microscope. Micro-positioning was carried out by a closed-loop three-axis piezo stage (Physik Instrumente) and AFM feedback control was carried out by a lock-in amplifier (HF2LI, Zurich Instruments). Photoluminescence of the NV centres was measured with an avalanche photodiode (Excelitas) and data were collected by a data acquisition card (PCIe-6353, National Instruments). Microwave pulses and sequences were created with a signal generator (Quicksyn FSW-0020, National Instruments) and modulated with an IQ mixer (Marki) and an arbitrary waveform generator (HDAWG, Zurich Instruments). NV centres were illuminated at <100 μW by a custom-designed 520 nm pulsed diode laser. Scanning NV tips were purchased from QZabre AG. Two tips were used throughout this study. Tip no. 1 (Figs. 2 and 3) had a stand-off distance of z = 95 ± 1 nm and tip no. 2 (Fig. 4) had a stand-off distance of z = 61 ± 2 nm (not including the 20 nm retract distance used while imaging), determined with a magnetic stripe41. Oscillation amplitudes were ~46 nm in Figs. 2b,d and 3e, ~92 nm in Fig. 2c and ~52 nm in Fig. 4a, determined via stroboscopic imaging12. A movable permanent neodymium magnet (supermagnet) below the sample served as the bias field source. Field alignment was possible by applying a fitting algorithm, which used a numerical model of the field produced by the magnet and NV resonance frequencies as a function of magnet position. The drift stability of the microscope was <30 nm per day and no drift correction techniques were employed.

Spin energy levels in a transverse magnetic field

With an off-axis field, the usual spin-state description (using ms) of the NV centre is not ideal, as the magnetic quantum number ms is no longer conserved. In this scenario the eigenstates, which we denote as \(\left\vert 0\right\rangle\), \(\left\vert -\right\rangle\) and \(\left\vert +\right\rangle\), are superpositions of the usual \(\left\vert {m}_{{{{\rm{s}}}}}\right\rangle\) states13,42. Spin-state mixing results in a reduced photoluminescence and optical contrast, which worsened the overall measurement sensitivity. This effect, however, is manageable for relatively weak (<12 mT) bias fields43. Bias fields <5 mT are also non-ideal, as 15N hyperfine coupling effectively misaligns small bias fields.

Alignment of the transverse magnetic field

The alignment algorithm of the bias field provided the ability to align with roughly 1° of uncertainty; however, we improved the alignment by additionally considering hyperfine interactions42. In an off-axis bias field (θB = 90°) the \(\left\vert 0\right\rangle\) state splits into two states that differ by Δ = 2ahfγeB/D, where ahf = 3.65 MHz is the off-axis 15N hyperfine coupling parameter, γe = 2π × 28 GHz T−1 is the gyromagnetic ratio of the electron, B is the applied magnetic field and D = 2.87 GHz is the zero-field splitting25. When deviating from θB = 90°, the on-axis component splits each of the \(\left\vert -\right\rangle\) and \(\left\vert +\right\rangle\) states through 15N hyperfine interaction resulting in eight total transitions (four for ω− and four for ω+). We swept the fitted polar angle a few degrees around 90° (with the same field magnitude) and tracked the ω− hyperfine resonances. The best alignment was achieved when only two resonances were visible and when the resonance frequency was the largest (on-axis contributions decreased the resonance frequency).

Influence of crystal strain

Internal strain in the diamond acts equivalently as a permanent d.c. E field via the piezocoupling coefficient44. For in-plane strains much weaker than electric signals (which is assumed during our analysis) the effect of strain is unimportant. For large in-plane strains (relative to the electric signal), it is still possible to carry out scanning gradiometry without a loss in sensitivity. In this case, the detection axis is controlled by the strain direction (Supplementary Section 1). The angular dependence changes from \(\cos \left(2{\varphi }_{B}+{\varphi }_{{E}_{{{{\rm{ac}}}}}}\right)\) for small strains to \(\cos ({\varphi }_{\xi }-{\varphi }_{{E}_{{{{\rm{ac}}}}}})\cos \left(2{\varphi }_{B}+{\varphi }_{\xi }\right)\), where φξ is the in-plane strain angle.

Laboratory and NV centre frames of reference

To translate between the laboratory frame and NV centre frame a representation of the NV centre’s crystallographic coordinate system is determined in the laboratory frame. The laboratory frame is defined by the vectors \(\hat{x}=[1,0,0]\), \(\hat{y}=[0,1,0]\) and \(\hat{z}=[0,0,1]\). The unit vector along the symmetry axis of the NV centre (\({\hat{z}}_{{{{\rm{NV}}}}}\)) pointing from the nitrogen atom towards the vacancy site, is chosen as the [111] crystallographic direction. The x unit vector (\({\hat{x}}_{{{{\rm{NV}}}}}\)) is taken to be orthogonal to \({\hat{z}}_{{{{\rm{NV}}}}}\) and pointing from the vacancy site towards one of the three nearest carbon atoms (\([11\overline{2}]\) for example, although there are three possible choices)13,45. Then \({\hat{y}}_{{{{\rm{NV}}}}}={\hat{z}}_{{{{\rm{NV}}}}}\times {\hat{x}}_{{{{\rm{NV}}}}}\) to preserve right-handedness. To translate between the crystallographic frame and the laboratory frame the bias field alignment algorithm and known (001) cut of the diamond tip is used. For example, with the NV centre angles of θNV = 55° and φNV = 0° (as shown in Fig. 1b), we get \({\hat{x}}_{{{{\rm{NV}}}}}=\frac{1}{\sqrt{3}}[1,0,-\sqrt{2}]\), \({\hat{y}}_{{{{\rm{NV}}}}}=[0,1,0]\) and \({\hat{z}}_{{{{\rm{NV}}}}}=\frac{1}{\sqrt{3}}[\sqrt{2},0,1]\). It is important to note that using a different NV centre reference frame definition can result in an incorrect computation of laboratory-frame electric fields (Supplementary Section 5).

Gradiometry technique

A complete description of scanning gradiometry, including calibration procedures, can be found in ref. 12. An ~2 μs laser pulse was used to polarize the NV centre into the ms = 0 state, which for small off-axis bias fields corresponds to the \(\left\vert 0\right\rangle\) state. Next, a microwave π/2 pulse is applied to create a superposition between \(\left\vert 0\right\rangle\) and one of \(\left\vert \pm \right\rangle\). The quantum phase ϕ accumulated between the two states during the coherent precession is \({\phi }_{\pm }=\int\nolimits_{0}^{\tau }g(t){{\Delta }}{\omega }_{\pm }(t){\mathrm{d}}t\), where g(t) is the modulation function46, \({{\Delta }}{\omega }_{\pm }(t)=\mp 2\uppi {k}_{\perp }{E}_{{{{\rm{ac}}}}}(t)\cos (2{\varphi }_{B}+{\varphi }_{{E}_{{{{\rm{ac}}}}}})\sin (2\uppi ft)\) is the detuning (see Supplementary Section 1) and τ is the evolution time. We used a four-phase cycling technique12,47 of the last π/2 pulse to measure ϕ±. The readout of the NV centre’s spin state was performed by another ~2 μs laser pulse, during which the photons emitted from the NV centre were collected across a ~600 ns window.

Samples

Lead zirconate titanate

The 50 nm thick Pb[Zr0.2Ti0.8]O3 film and the 10 nm thick SrRuO3 electrode were grown on (001)-oriented SrTiO3 (CrysTec) using pulsed-layer deposition with a KrF excimer laser at 248 nm (LPXpro, Coherent). SrRuO3 was grown at a substrate temperature of 700 °C with an O2 partial pressure of 0.1 mbar and a laser fluence of 0.95 J cm−2 at 4 Hz. Pb[Zr0.2Ti0.8]O3 was grown at 550 °C at 0.12 mbar O2 partial pressure and a laser fluence of 1.2 J cm−2 at 4 Hz. The film was subsequently cooled to room temperature under growth pressure. Layer thicknesses were measured using X-ray reflectivity with a four-cycle thin-film diffractometer (PANalytical X’Pert3 MRD, CuKα1). Topography and PFM experiments were performed on a Bruker Multimode 8 atomic force microscope using Pt-coated Si tips (MikroMasch, k = 5.4 N m−1).

Hexagonal yttrium manganite

The YMnO3 bulk crystal was grown by the floating-zone technique, pre-oriented using Laue diffraction and cut perpendicular to the crystal z axis with a diamond saw. The sample was flattened by lapping with Al2O3 powder in water solution (9 μm particle size). Subsequently, the sample was chemomechanically polished using a colloidal silica slurry. To generate domains, the sample was pre-annealed and cooled through the Curie temperature TC in an O2 atmosphere48.

Electric field vector reconstruction

To reconstruct the E field vector (Fig. 2f and Extended Data Fig. 4) the two images recorded with a 45° difference in φB, denoted as \({I}_{1}(x,y)={E}_{{{{\rm{ac}}}}}\cos (2{\varphi }_{B}+{\varphi }_{{E}_{{{{\rm{ac}}}}}})\) and \({I}_{2}(x,y)={E}_{{{{\rm{ac}}}}}\sin (2{\varphi }_{B}+{\varphi }_{{E}_{{{{\rm{ac}}}}}})\), are combined to yield the magnitude (Eac(x, y)) and angle (\({\varphi }_{{E}_{{{{\rm{ac}}}}}}(x,y)\)) of the measured electric field signal

$${E}_{{{{\rm{ac}}}}} (x,y)=\sqrt{{({I}_{1}(x,y))}^{2}+{({I}_{2}(x,y))}^{2}},$$
(2)
$${\varphi }_{{E}_{{{{\rm{ac}}}}}}(x,y)=\arctan \left(\frac{{I}_{2}(x,y)}{{I}_{1}(x,y)}\right)-2{\varphi }_{B}\,(\,{{\mbox{mod}}}\,\,2\uppi ).$$
(3)

Since there are three possible choices for \({\hat{x}}_{{{{\rm{NV}}}}}\), owing to the C3ν symmetry of the NV centre, φB and \({\varphi }_{{E}_{{{{\rm{ac}}}}}}(x,y)\) are only known up to a multiple of 2π/3. This propagates to the E field gradients along the x and y directions of the NV centre, which are calculated as

$${x}_{{{{\rm{osc}}}}}{\partial }_{r}{E}_{{x}_{{{{\rm{NV}}}}}}(x,y)= {E}_{{{{\rm{ac}}}}} (x,y)\cos ({\varphi }_{{E}_{{{{\rm{ac}}}}}}(x,y)),$$
(4)
$${x}_{{{{\rm{osc}}}}}{\partial }_{r}{E}_{{y}_{{{{\rm{NV}}}}}}(x,y)= {E}_{{{{\rm{ac}}}}} (x,y)\sin ({\varphi }_{{E}_{{{{\rm{ac}}}}}}(x,y)),$$
(5)

where xosc is the tip oscillation amplitude and ∂r is the directional derivative along the unit vector \(\hat{r}=\hat{x}\cos \alpha +\hat{y}\sin \alpha\) in the laboratory frame and α is the in-plane oscillation angle. With either of these images it is possible to reconstruct the laboratory components of the E field gradients in Fourier space. For example, with \({\hat{x}}_{{{{\rm{NV}}}}}\) the laboratory-frame E field gradient vector is

$${{{\mathcal{F}}}}\{{\partial }_{r}{{{{\bf{E}}}}}^{{{{\rm{lab}}}}}\}={{{\mathcal{F}}}}\{{\partial }_{r}{E}_{{x}_{{{{\rm{NV}}}}}}\}\frac{{{{\bf{K}}}}}{{\hat{x}}_{{{{\rm{NV}}}}}\cdot {{{\bf{K}}}}},$$
(6)

where \({{{\mathcal{F}}}}\) is the Fourier transform, K = [ikx, iky, K] and \(K=\sqrt{{k}_{x}^{2}+{k}_{y}^{2}}\). The laboratory-frame vector components can be determined independent of the three possible \({\hat{x}}_{{{{\rm{NV}}}}}\) and \({\hat{y}}_{{{{\rm{NV}}}}}\) choices because the term in the denominator removes its influence. The last step is integration in Fourier space with a wavevector-dependent window function that cuts off high-frequency terms (using a Hann window filter)47, and a line filter that removes the amplified noise perpendicular to the direction of integration12. The E field vector is computed with

$${{{\mathcal{F}}}}\{{{{{\bf{E}}}}}^{{{{\rm{lab}}}}}\}=\frac{{{{\mathcal{F}}}}\{{\partial }_{r}{{{{\bf{E}}}}}^{{{{\rm{lab}}}}}\}W(\lambda ,\alpha )}{-i{x}_{{{{\rm{osc}}}}}{k}_{r}},$$
(7)

where \({k}_{r}={k}_{x}\cos (\alpha )+{k}_{y}\sin (\alpha )\) and W(λ, α) is the window function. We set the cut-off wavelength to the stand-off distance (λ = z, which produces a cut-off wavevector of k = 2π/λ) and the oscillation angle to match the x direction (α = 0°). We applied this procedure on both the \({\partial }_{x}{E}_{{x}_{{{{\rm{NV}}}}}}\) and \({\partial }_{x}{E}_{{y}_{{{{\rm{NV}}}}}}\) images and average the results.

Surface charge density reconstruction

To reverse propagate our E field measurements into a surface charge density (Fig. 2g), we first treat Coulomb’s law for a two-dimensional surface charge density σ(x,y) as a convolution integral in Fourier space49. With the transfer function G(K, z) = e−Kz/(2ϵ0K), where ϵ0 is the vacuum permittivity, the E field components from a single surface charge density become

$${{{\mathcal{F}}}}\{{{{{\bf{E}}}}}^{{{{\rm{lab}}}}}\}={{{\mathcal{F}}}}\{\sigma \}G(K,z){{{\bf{K}}}}.$$
(8)

For two surface charge densities of opposite polarity separated by a distance t, the Fourier transformed E field is \({{{\mathcal{F}}}}\{{{{{\bf{E}}}}}^{{{{\rm{total}}}}}\}={{{\mathcal{F}}}}\{{{{{\bf{E}}}}}^{{{{\rm{lab}}}}}\}(1-{\mathrm{e}}^{-Kt})\), where the −e−Kt term comes from the bottom surface. The surface charge density can be computed for each of the three laboratory-frame vector components, and from our NV electrometry measurements the surface charge density is

$${{{\mathcal{F}}}}\{\sigma \}=\frac{2{\epsilon }_{0}K}{3{x}_{{{{\rm{osc}}}}}{k}_{r}}\frac{{\mathrm{e}}^{Kz}}{1-{\mathrm{e}}^{-Kt}}\left(\frac{{{{\mathcal{F}}}}\{{\partial }_{r}{E}_{x}^{{{{\rm{lab}}}}}\}}{{k}_{x}}+\frac{{{{\mathcal{F}}}}\{{\partial }_{r}{E}_{y}^{{{{\rm{lab}}}}}\}}{{k}_{y}}+\frac{i{{{\mathcal{F}}}}\{{\partial }_{r}{E}_{z}^{{{{\rm{lab}}}}}\}}{K}\right),$$
(9)

where the three vector components have been averaged. We apply an additional window function, W(λ1, λ2, α), to equation (9) that cuts off both low-frequency (λ1 = 30z) and high-frequency (λ2 = z) components, and a line filter (α = 0°) to remove amplified noise from the deconvolution process.

Electric field gradiometry line scan fitting

Fitting the line scans in Fig. 3e was accomplished by first determining a simplified form for different surface charge models (Supplementary Section 2). For a monopole domain wall located at x = xi and propagating along y, the equations

$$\begin{array}{l}{\partial }_{x}{E}_{x}^{{{{\rm{lab}}}}}=\frac{-\sigma }{\uppi {\epsilon }_{0}}\frac{x-{x}_{i}}{{(x-{x}_{i})}^{2}+{z}^{2}},\\ {\partial }_{x}{E}_{z}^{{{{\rm{lab}}}}}=\frac{-\sigma }{\uppi {\epsilon }_{0}}\frac{z}{{(x-{x}_{i})}^{2}+{z}^{2}},\end{array}$$
(10)

are used. For the equivalent dipole domain wall the equations

$$\begin{array}{l}{\partial }_{x}{E}_{x}^{{{{\rm{lab}}}}}=\frac{\sigma d}{\uppi {\epsilon }_{0}}\frac{2(x-{x}_{i})z}{{({(x-{x}_{i})}^{2}+{z}^{2})}^{2}},\\ {\partial }_{x}{E}_{z}^{{{{\rm{lab}}}}}=\frac{\sigma d}{\uppi {\epsilon }_{0}}\frac{{z}^{2}-{(x-{x}_{i})}^{2}}{{({(x-{x}_{i})}^{2}+{z}^{2})}^{2}},\end{array}$$
(11)

are used. In both types of domain walls \({\partial }_{x}{E}_{y}^{{{{\rm{lab}}}}}=0\). Here, the stand-off z, surface charge density σ (or surface dipole density σd), domain wall locations (x1 and x2) and sample thickness t (for the bottom layer) are used to create the different surface charge models. Next, the gradient components are projected onto the NV centre’s x and y unit vectors using \({\partial }_{x}{E}_{{x}_{{{{\rm{NV}}}}}}={\partial }_{x}{{{{\bf{E}}}}}^{{{{\rm{lab}}}}}\cdot {\hat{x}}_{{{{\rm{NV}}}}}\) and \({\partial }_{x}{E}_{{y}_{{{{\rm{NV}}}}}}={\partial }_{x}{{{{\bf{E}}}}}^{{{{\rm{lab}}}}}\cdot {\hat{y}}_{{{{\rm{NV}}}}}\). The unit vectors depend on the polar and azimuthal angles θNV and φNV. Then, the E field angle is computed as \({\varphi }_{{\mathrm{E}}_{{{{\rm{ac}}}}}}=\arctan ({\partial }_{x}{E}_{{y}_{{{{\rm{NV}}}}}}/{\partial }_{x}{E}_{{x}_{{{{\rm{NV}}}}}})\). Finally, the measured signal is modelled by \({E}_{{{{\rm{meas}}}}}={x}_{{{{\rm{osc}}}}}\sqrt{{\left({\partial }_{x}{E}_{{x}_{{{{\rm{NV}}}}}}\right)}^{2}+{\left({\partial }_{x}{E}_{{y}_{{{{\rm{NV}}}}}}\right)}^{2}}\cos (2{\varphi }_{B}+{\varphi }_{{E}_{{{{\rm{ac}}}}}})\). During the fitting process, the NV angles, magnetic bias field angle, sample thickness, stand-off distance and oscillation amplitude are kept constant, having been measured or determined previously. The surface charge density (or surface dipole density) and domain wall locations are fitted and the best model is determined by the shape, quality and polarity of the fit.

Estimation of sensitivity

We estimated the sensitivity with two methods, first by error propagating the measurement counts used in the quantum phase computation, and second by taking line-by-line differences from two consecutive line scans and computing the standard deviation of the resulting trace47. As shown in Supplementary Section 3, our best sensitivities were achieved by using multipulse sequences with quantum phase accumulation across multiple oscillation periods12. The error propagated sensitivity was 0.24 kV cm−1 Hz−1/2 and the line-by-line sensitivity was 0.29 kV cm−1 Hz−1/2.