Imaging ferroelectric domains with a single-spin scanning quantum sensor

The ability to sensitively image electric fields is important for understanding many nanoelectronic phenomena, including charge accumulation at surfaces1 and interfaces2 and field distributions in active electronic devices3. A particularly exciting application is the visualization of domain patterns in ferroelectric and nanoferroic materials4,5, owing to their potential in computing and data storage6–8. Here, we use a scanning nitrogen-vacancy (NV) microscope, well known for its use in magnetometry9, to image domain patterns in piezoelectric (Pb[Zr0.2Ti0.8]O3) and improper ferroelectric (YMnO3) materials through their electric fields. Electric field detection is enabled by measuring the Stark shift of the NV spin10,11 using a gradiometric detection scheme12. Analysis of the electric field maps allows us to discriminate between different types of surface charge distributions, as well as to reconstruct maps of the three-dimensional electric field vector and charge density. The ability to measure both stray electric and magnetic fields9,13 under ambient conditions opens opportunities for the study of multiferroic and multifunctional materials and devices8,14.

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 microscopy 15 and piezoresponse force microscopy (PFM) 4 , along with the related techniques of Kelvin probe force microscopy 16 , low-energy electron microscopy 5 and emerging scanning quantum technologies 1,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 diamond 9,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 tip 12 to overcome the challenges of static screening 19 and low coupling to electric fields 10,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 diamond 11 . We illustrate the impact of our approach by imaging patterned domains in application-relevant 6 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 (z NV axis). This leads to an in-plane electric field coupling that is approximately 50× larger than the out-of-plane coupling 11,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 electrodes 11,20,21 , charged scanning probes 22,23 , surface band bending 24 and intrinsic dopant charges 25 . Letter https://doi.org/10.1038/s41567-022-01921-4 This issue has hindered previous attempts at implementing a scanning NV electrometer 19,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 direction 12 . In addition to alleviating static screening, the a.c. detection also improves sensitivity by at least an order of magnitude 11,12,18 . For the spin-echo pulse sequence shown in Fig. 1d, the field-induced coherent phase accumulation of the NV spin is 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[Zr 0.2 Ti 0.8 ]O 3 , PZT) grown on top of a SrRuO 3 -buffered (001)-oriented SrTiO 3 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. 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 signatures 12 . 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 state 27 . 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 PZT 28 . 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 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).
To enable electric field (E) detection, we accurately orient an external bias field of 5−12 mT transverse to z NV (θ B = 90°) with ~0.5° of uncertainty (Methods). In this bias field configuration, the spin transition frequencies ω ± are linearly sensitive to electric fields, while correcting for magnetic fields up to second order 11,13 . Here, ω (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 constant 21 , E ⊥ is the magnitude of the in-plane E field vector and φ B = arctan(B y NV /B x NV ) and φ E = arctan(E y NV /E x 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). 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 lab z 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 surfaces 29 , even once the materials are buried in a device architecture 30 .
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 adsorbates 32 or the formation of a charged, off-stoichiometric surface layer 33 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 screening 34 .
We further illustrate scanning electrometry by imaging the natural domain pattern of an improper ferroelectric, the hexagonal manganite YMnO 3 . Hexagonal manganites are exciting benchmark materials because their surface polarization (P ~ 5.5 μC cm −2 for YMnO 3 ) is over an order of magnitude smaller than our PZT sample and typical for ferroelectric and multifunctional materials 8,14 . In addition, hexagonal manganites are type I multiferroics and become antiferromagnetically ordered below T N ~ 100 K (refs. 14,35 ). While the antiferromagnetic domain pattern has been imaged with scanning NV magnetometry at cryogenic temperatures 26 , here, we focus on the ferroelectric domain pattern accessible under ambient conditions. Figure 4a shows an electrometry map recorded above a polished, bulk YMnO 3 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).
In summary, we demonstrate nanoscale electric field imaging using a scanning NV microscope. Crucial to our experiment is the use of gradiometric detection 12 , 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 bulk 11 .
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 techniques 4,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 temperatures 37 . 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 materials 35,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 magnetization 39,40 .
Note added in proof: While finalizing this manuscript, we became aware of related work 34 describing scanning NV electrometry using a similar gradiometry technique to image electric fields from electrodes.

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/s41567-022-01921-4. x NV x NV

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 stripe 41 .
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 m s ) of the NV centre is not ideal, as the magnetic quantum number m s is no longer conserved. In this scenario the eigenstates, which we denote as |0⟩, |−⟩ and |+⟩, are superpositions of the usual |m s ⟩ states 13,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 fields 43 . Bias fields <5 mT are also non-ideal, as 15 N 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 interactions 42 . In an off-axis bias field (θ B = 90°) the |0⟩ state splits into two states that differ by Δ = 2a hf γ e B/D, where a hf = 3.65 MHz is the off-axis 15 N 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 splitting 25 . When deviating from θ B = 90°, the on-axis component splits each of the |−⟩ and |+⟩ states through 15 N 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 coefficient 44 . 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 (2φ B + φ E ac ) for small strains to 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 x = [1, 0, 0], ŷ = [0, 1, 0] and ẑ = [0, 0, 1]. The unit vector along the symmetry axis of the NV centre (ẑ NV ) pointing from the nitrogen atom towards the vacancy site, is chosen as the [111] crystallographic direction. The x unit vector (x NV ) is taken to be orthogonal to ẑ NV and pointing from the vacancy site towards one of the three nearest carbon atoms ([112] for example, although there are three possible choices) 13,45 . Then ŷ NV =ẑ NV ×x 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 x

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 m s = 0 state, which for small off-axis bias fields corresponds to the |0⟩ state. Next, a microwave π/2 pulse is applied to create a superposition between |0⟩ and one of |±⟩. The quantum phase ϕ accumulated between the two states during the coherent precession is ϕ ± = ∫ τ 0 g(t)Δω ± (t)dt, where g(t) is the modulation function 46 , Δω ± (t) = ∓2 k ⟂ E ac (t) cos(2φ B + φ E ac ) sin(2 ft) is the detuning (see Supplementary Section 1) and τ is the evolution time. We used a four-phase cycling technique 12,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. Hexagonal yttrium manganite. The YMnO 3 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 Al 2 O 3 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 T C in an O 2 atmosphere 48 .

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 ac cos(2φ B + φ E ac ) and I 2 (x, y) = E ac sin(2φ B + φ E ac ), are combined to yield the magnitude (E ac (x, y)) and angle (φ E ac (x, y)) of the measured electric field signal Letter https://doi.org/10.1038/s41567-022-01921-4 φ E ac (x, y) = arctan ( I 2 (x, y) Since there are three possible choices for x NV , owing to the C 3ν symmetry of the NV centre, φ B and φ E 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 osc ∂ r E y NV (x, y) = E ac (x, y) sin(φ E ac (x, y)), (5) where x osc is the tip oscillation amplitude and ∂ r is the directional derivative along the unit vector r =x cos α +ŷ sin α 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 x NV the laboratory-frame E field gradient vector is where ℱ is the Fourier transform, K = [ik x , ik y , K] and K = √ k 2 x + k 2 y . The laboratory-frame vector components can be determined independent of the three possible x NV and ŷ 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 integration 12 . The E field vector is computed with where k r = k x cos(α) + k y sin(α) 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 ∂ x E x NV and ∂ x E y 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 space 49 .
With the transfer function G(K, z) = e −Kz /(2ϵ 0 K), where ϵ 0 is the vacuum permittivity, the E field components from a single surface charge density become For two surface charge densities of opposite polarity separated by a distance t, the Fourier transformed E field is ℱ{E total } = ℱ{E lab }(1 − 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 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 = x i and propagating along y, the equations are used. For the equivalent dipole domain wall the equations are used. In both types of domain walls ∂ x E lab y = 0. Here, the stand-off z, surface charge density σ (or surface dipole density σd), domain wall locations (x 1 and x 2 ) 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 ∂ x E x NV = ∂ x E lab ⋅x NV and ∂ x E y NV = ∂ x E lab ⋅ŷ NV . The unit vectors depend on the polar and azimuthal angles θ NV and φ NV . Then, the E field angle is computed as φ E ac = arctan(∂ x E y NV /∂ x E x NV ). Finally, the measured signal is modelled by E meas = x osc√ (∂ x E x NV ) 2 + (∂ x E y NV ) 2 cos(2φ B + φ E 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 trace 47 . As shown in Supplementary Section 3, our best sensitivities were achieved by using multipulse sequences with quantum phase accumulation across multiple oscillation periods 12 . 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 .

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.