Spatial mapping of band bending in semiconductor devices using in-situ quantum sensors

Spatial mapping of band bending in semiconductor devices using in-situ quantum sensors

D. A. Broadway School of Physics, University of Melbourne, Parkville, VIC 3010, Australia Centre for Quantum Computation and Communication Technology, School of Physics, University of Melbourne, Parkville, VIC 3010, Australia    N. Dontschuk School of Physics, University of Melbourne, Parkville, VIC 3010, Australia Centre for Quantum Computation and Communication Technology, School of Physics, University of Melbourne, Parkville, VIC 3010, Australia    A. Tsai School of Physics, University of Melbourne, Parkville, VIC 3010, Australia    S. E. Lillie School of Physics, University of Melbourne, Parkville, VIC 3010, Australia Centre for Quantum Computation and Communication Technology, School of Physics, University of Melbourne, Parkville, VIC 3010, Australia    C. T.-K. Lew School of Physics, University of Melbourne, Parkville, VIC 3010, Australia Centre for Quantum Computation and Communication Technology, School of Physics, University of Melbourne, Parkville, VIC 3010, Australia    J. C. McCallum School of Physics, University of Melbourne, Parkville, VIC 3010, Australia    B. C. Johnson School of Physics, University of Melbourne, Parkville, VIC 3010, Australia Centre for Quantum Computation and Communication Technology, School of Physics, University of Melbourne, Parkville, VIC 3010, Australia    M. W. Doherty Laser Physics Centre, Research School of Physics and Engineering, Australian National University, Canberra, ACT 2601, Australia    A. Stacey Centre for Quantum Computation and Communication Technology, School of Physics, University of Melbourne, Parkville, VIC 3010, Australia Melbourne Centre for Nanofabrication, Clayton, VIC 3168, Australia    L. C. L. Hollenberg School of Physics, University of Melbourne, Parkville, VIC 3010, Australia Centre for Quantum Computation and Communication Technology, School of Physics, University of Melbourne, Parkville, VIC 3010, Australia    J.-P. Tetienne School of Physics, University of Melbourne, Parkville, VIC 3010, Australia
thanks: These authors contributed equally to this work.thanks: These authors contributed equally to this work.

Band bending is a central concept in solid-state physics that arises from local variations in charge distribution especially near semiconductor interfaces and surfaces Zhang2012 (); Kotadiya2018 (); Simon2010 (). Its precision measurement is vital in a variety of contexts from the optimisation of field effect transistors Stathis2006 (); Zhang2006 (); Kaczr2018 () to the engineering of qubit devices with enhanced stability and coherence Weber2010 (); Kaviani2014 (); Usman2016 (). Existing methods are surface sensitive and are unable to probe band bending at depth from surface or bulk charges related to crystal defects Zhang2012 (); Kronik1999 (); Ishii2004 (); Butler2014 (). Here we propose an in-situ method for probing band bending in a semiconductor device by imaging an array of atomic-sized quantum sensing defects to report on the local electric field. We implement the concept using the nitrogen-vacancy centre in diamond Doherty2013 (); Dolde2011 (), and map the electric field at different depths under various surface terminations. We then fabricate a two-terminal device based on the conductive two-dimensional hole gas formed at a hydrogen-terminated diamond surface Strobel2004 (), and observe an unexpected spatial modulation of the electric field attributed to a complex interplay between charge injection and photo-ionisation effects. Our method opens the way to three-dimensional mapping of band bending in diamond and other semiconductors hosting suitable quantum sensors, combined with simultaneous imaging of charge transport in complex operating devices Tetienne2017 ().

The emergence of semiconductor-based quantum sensing technologies in the last decade has opened new opportunities in a range of disciplines across physics, materials science and biology Degen2017 (). While most existing applications involve sensors that are external to the target sample to be measured Rondin2014 (); Casola2018 (), in-situ quantum sensors can also be an extremely valuable resource to study the sample itself by enabling three-dimensional (3D) mapping Iwasaki2017 (). For semiconductor materials this is especially advantageous as it allows information to be gained on the interactions between surface and bulk defects, which play an important role in semiconductor electronics and in quantum technologies. Here we propose and demonstrate such an application, where in-situ quantum sensors are used to map the electric field near a semiconductor surface, , which is related to band bending via the standard relation


where is the energy of electrons at the valence band maximum relative to the Fermi level ( is the electron charge). Specifically, we employ the nitrogen-vacancy (NV) centre in diamond Doherty2013 (), a well established atomic-sized quantum sensing system that was recently shown to be a sensitive electrometer Dolde2011 (); Dolde2014 (), although the concept could be applied to other semiconductor systems such as silicon Zhang2018 () or silicon carbide Falk2014 (). NV centres can be positioned with a resolution below 1 nm in one dimension and 20 nm in the other two Ohno2012 (); Lesik2016 (), making it an ideal candidate for 3D mapping of built-in electric fields.

Figure 1: Mapping band bending with in-situ quantum sensors. a, Principle of the experiment, where NV sensors (green dots) probe the electric field associated with surface band bending, here visualised as a distribution of space charge density. b, Calculated electric field profile for a typical (001)-oriented, oxygen-terminated diamond surface, modelled as a layer of surface acceptor defects with a density of states centred at  eV above the valence band maximum Stacey2018 () and a surface density  nm (see details in SI). The implanted substitutional nitrogen and NV defects are taken to be uniformly distributed over the range  nm (i.e.  nm), with a total areal density of  nitrogen/nm. The green line is the NV population at equilibrium, and the green (red) shading represents the region where the NV (NV) charge state is dominant. Inset: corresponding band diagram near the surface, where is the conduction band minimum, the valence band maximum, the Fermi level, and NV represents the charge transition level of the NV centre, i.e. NV is the stable charge state when this level is below . c, Diagram of the experimental set-up showing the diamond sample mounted on a glass slide patterned with gold to allow microwave (MW) injection and interfacing with electrical devices, illumination with a green laser and imaging of the NV red photoluminescence (PL) with a camera. d, Example ODMR spectrum recorded for an ensemble of near-surface NV centres in an oxygen-terminated diamond under a magnetic field chosen to be perpendicular (within less than ) to a given NV family (here NV). Each resonance is labelled according to the corresponding NV orientation, defined in e. The solid line is a multiple-Lorentzian fit. e, The four possible tetrahedral orientations of the NV bond with respect to the sample reference frame .

The principle of the experiment is depicted in Fig. 1a. At the diamond surface, the bands bend to neutralise any surface charge due to ionised adsorbates or surface defect states, resulting in an electric field perpendicular to the surface with a magnitude . To probe this electric field, nitrogen ions were implanted to form NV centres, following a spatial distribution that can be approximated as uniform over the depth range , where is the (tunable) mean implantation depth Lehtinen2016 (). To estimate , we first consider the case of commonly used oxygen-terminated diamond. It was recently found that such samples typically host surface defects that introduce an acceptor level into the band-gap, with densities () as high as 1 nm Stacey2018 (). An example of a calculated electric field profile for this scenario (with parameters representative of our implanted samples) is plotted in Fig. 1b, predicting a maximum value at the surface of  MV/cm and a characteristic decay length of  nm. We note that is positive (i.e., the electric field points towards the surface) which corresponds to the bands bending upward (inset in Fig. 1b), as expected from a positive space charge density near the surface (see Fig. 1a). As a consequence, only NVs deeper than a certain threshold (here  nm for  nm) exist in the negatively charged state (NV) usable for sensing (Fig. 1b). The expectation value for an electric field measurement, i.e. averaged over the NV distribution, is  kV/cm, well in the range of sensitivity of the NV centre Dolde2011 ().

To measure this electric field, we performed optically detected magnetic resonance (ODMR) spectroscopy on the NV centres, using the experimental set-up depicted in Fig. 1c. An example ODMR spectrum obtained from a near-surface NV ensemble ( nm) is shown in Fig. 1d. A small magnetic field (of magnitude  mT) was applied to split the eight otherwise degenerate electron spin resonances corresponding to the four possible NV defect orientations relative to the diamond crystal (Fig. 1e), and carefully oriented so as to maximise the sensitivity to electric fields Dolde2011 () (see SI). The spectrum was fit to extract the eight resonance frequencies, which are then compared to the standard NV spin Hamiltonian including the Zeeman and Stark effects Dolde2011 (); Doherty2012 (), allowing us to infer the full vector magnetic and electric fields simultaneously (except for an overall sign ambiguity, i.e. and yield the same ODMR spectrum). Importantly, we found that the measured frequencies could be satisfactorily fit only when accounting for the Stark effect, providing clear evidence of the presence of a non-vanishing electric field (see fit error analysis in SI). For the data shown in Fig. 1d, we obtain  kV/cm (where we fixed to reduce the uncertainty), reasonably close to our estimate.

Figure 2: Electric field versus implantation depth and surface termination. a, Electric field, , as a function the mean implantation depth, , for O-terminated diamond. Solid lines: result of the band bending model described in Fig. 1b with  nm (lower curve) and  nm (upper), with the shading representing intermediate values. Inset: Difference between the electric field measured for H- and O-terminated diamond. Solid lines: model using a fixed density of charged surface adsorbates,  nm, and with  nm (upper curve) and  nm (lower). b, map of an H-terminated channel on an O-terminated background ( nm). c, PL reduction of the H region relative to the O region, as a function of . Solid lines: model using the same parameters as in inset of a. Inset: PL image of an H-terminated channel ( nm). d, vs calculated for  nm. The dashed lines and data points indicate the measured values for a comparable sample following various surface treatments, performed in order: oxygen plasma (as used to form the O termination in a-c), acid cleaning, oxygen burning and piranha treatment (see details in SI). Vertical error bars: where is the standard deviation.

To illustrate the 3D mapping capability, we formed NV centres at different depths in distinct diamonds and measured the electric field as explained above. We found that decreased from  kV/cm at  nm to  kV/cm at  nm (Fig. 2a), in good agreement with our modelling for in the range  nm. We next studied the effect of surface termination on the electric field by forming a hydrogen-terminated (H) channel on an otherwise oxygen-terminated (O) diamond. An example of the resulting map is shown in Fig. 2b, revealing an increase from  kV/cm in the O region to  kV/cm in the H region. This is expected because H-terminated diamond is known to have a lower electron affinity, leading to efficient charge transfer from the diamond material onto acceptor species adsorbed on the surface in ambient air Strobel2004 (). This leads to increased band bending and hence the increased electric field, beyond the threshold required to form a conductive two-dimensional hole gas (2DHG) near the surface Strobel2004 (); Pakes2014 (), which is imaged in an electrical device described below. By fitting our model to the measured increase in caused by the H termination (inset in Fig. 2a), one can infer the density of charged surface adsorbates (acceptors),  nm, in good agreement with the value derived from surface resistivity measurements (see SI).

A consequence of the increased band bending is a decrease in the number of NV centres, hence a decrease in detected photoluminescence (PL), since the NVs closest to the surface become charge neutral Hauf2011 (). This is illustrated in Fig. 2c, which shows the PL reduction as a function of , with an example PL image of a H channel shown in inset. This motivates the need to minimise band bending via surface engineering for quantum sensing applications, where the NV to surface distance is critical Rondin2014 (). As a step towards this goal, we applied various surface treatments in an attempt to reduce the density of surface defects. Starting from a diamond initially O-terminated with an oxygen plasma as previously ( kV/cm), we were able to reduce the electric field to  kV/cm through a combination of wet and dry treatments, which corresponds to a reduction of by nearly a factor 3 according to our theory (Fig. 2d), and a reduction in the mean NV depth from  nm to  nm. These trends are broadly consistent with direct measurements of reported recently Stacey2018 (). We note that another avenue to reduce is by etching the diamond, as shown in the SI. These results illustrate the value of in-situ electric field measurements, which provide new insights into semiconductor surfaces.

Figure 3: Electric field in a two-terminal device. a, Photograph of the diamond with the fabricated devices. b, Optical bright-field image of a device showing the TiC/Pt/Au contacts in dark, defined as source (S) and drain (D). The optically-transparent 2DHG channel is indicated by orange dotted lines. The black dashed box defines the area imaged in c-h. c,d map under an applied voltage  V (c) and  V (d), corresponding to a current through the channel of A. e, Line cut of the PL along the channel, indicative of the laser spot profile (symbolised by the shading under the curve). f-h, Same as c-e but with the laser spot offset by m to the right. i, Diagram illustrating the different recombination and trapping processes for electrons under laser illumination and bias voltage. The TiC contact extends  nm deep into the diamond (see SI), with the NV centres implanted at a mean depth  nm. j, Line cuts of taken along the channel (see dashed line in f) for different voltages , with the laser spot centred as in f-h. k, Electric field measured at the position indicated by the black dashed line in j, plotted as a function of . In j,k, the zero-voltage value was subtracted to show only the increase caused by the applied voltage, .

We now demonstrate mapping of the electric field in an active electrical device consisting of a driven conductive 2DHG channel formed with an H-terminated diamond surface. Two-terminal devices were fabricated where TiC/Pt/Au contacts (source and drain) are connected by a H-terminated channel 100 m in length and 20 m in width (Figs. 3a,b). Unexpectedly, upon applying a voltage  V an increase in was observed (by up to a factor 2) in a well-defined region of the channel extending over m from the drain (Fig. 3c). Upon inversion of the voltage ( V), the feature moved to the other contact (Fig. 3d) to remain at the hole drain. In addition, we observed an influence of the position of the laser illumination spot used for the measurements. In Figs. 3c,d, the laser spot was centred relative to the device as shown by the PL profile in Fig. 3e. When the laser was moved by m towards the right-hand contact (Figs. 3f-h), the region of increased electric field appeared to extend further away (m) from the left-hand contact under positive voltage (Fig. 3f), and disappeared completely under negative voltage (Fig. 3g).

These observations are qualitatively interpreted as a combination of two competing effects, illustrated in Fig. 3i: the injection of electrons from the drain contact into charge traps in the diamond bulk, increasing the electric field seen by the NV centres, and laser-induced ionisation of these charge traps, allowing the electrons to escape via the conduction band and returning the electric field to its zero-voltage value. This interpretation is corroborated by the negative capacitance measured for these devices (see SI). Line cuts taken at different voltages (Fig. 3j) reveal that the electric field decreases steadily from the contact (with a maximum value that increases monotonically with voltage, see Fig. 3k) before dropping off abruptly at a specific position independent of voltage (but dependent on laser position), even though the laser intensity increases approximately linearly with position (Fig. 3h). This suggests the existence of a strong non-linearity in the ionisation process as a function of laser intensity, possibly due to a change in the ionisation energy of the charge traps caused by ionisation of a second species of defects Dhomkar2018 (). These experiments illustrate how previously un-observable lateral changes in electric field, resulting from a complex contact/diamond junction, can be directly imaged and correlated with the electrical properties of the device (here a negative capacitance).

Looking forward, our method could be used to validate models of band bending used in semiconductor electronics, or to facilitate the optimisation of surfaces for semiconductor-based quantum technologies, including the improvement of the charge stability and quantum coherence of near-surface NV centres in diamond, which remains an outstanding problem to date DeOliveira2017 (); Stacey2018 (). Another exciting prospect is to combine electric field mapping with other quantum sensing modalities, e.g. current flow mapping Tetienne2017 () or noise spectroscopy Casola2018 (); Kolkowitz2015 (), which would allow simultaneous monitoring of charge transport and band bending, as in a three-terminal field-effect transistor device based on the native 2DHG formed on diamond or on other 2D electronic systems such as graphene.

Author contributions. NV measurements and analysis were performed by D.A.B and J.-P.T, with inputs from M.W.D. The devices were fabricated by N.D and D.A.B, H-terminated by A.T and A.S, and electrically characterised by C.T-K.L and B.C.J. The band bending model was developed by N.D with inputs from D.A.B, J.-P.T, A.S and L.C.L.H. All authors contributed to interpreting the data and writing the manuscript.

Acknowledgements. We thank Michael Barson, David Simpson and Liam Hall for useful discussions. We acknowledge support from the Australian Research Council (grants CE110001027, DE170100129, FL130100119, DP170102735). J.-P.T acknowledges support from the University of Melbourne through an Early Career Researcher Grant. D.A.B, A.T, S.E.L and C.T.-K.L are supported by an Australian Government Research Training Program Scholarship. This work was performed in part at the Melbourne Centre for Nanofabrication (MCN) in the Victorian Node of the Australian National Fabrication Facility (ANFF).

I Experimental methods

i.1 Diamond samples

The NV-diamond samples used in these experiments were made from 4 mm 4 mm 50 m electronic-grade ([N]  ppb) single-crystal diamond plates with {110} edges and a (001) top facet, purchased from Delaware Diamond Knives, subsequently laser cut into smaller 2 mm 2 mm 50 m plates and acid cleaned (15 minutes in a boiling mixture of sulphuric acid and sodium nitrate). Each plate was then implanted with N ions (InnovIon) at a given energy (ranging from 4 to 20 keV, see Table S1), a dose of ions/cm, and with a tilt angle of 7. Following implantation, the diamonds were annealed in a vacuum of  Torr to form the NV centres, using the following sequence Tetienne2018 (): 6h at 400C, 2h ramp to 800C, 6h at 800C, 2h ramp to 1100C, 2h at 1100C, 2h ramp to room temperature. To remove the graphitic layer formed during the annealing at the elevated temperatures, the samples were acid cleaned (15 minutes in a boiling mixture of sulphuric acid and sodium nitrate) before any further surface treatment or fabrication step.

Diamond Implantation energy Figures
#1 6 keV 21 nm 1d, 2a, 2b, 2c
#2 10 keV 35 nm 2a, 2c, 2d, 3, S5-10
#3 4 keV 14 nm 2a, 2c & inset
#4 14 keV 49 nm 2a, 2c
#5 20 keV 70 nm 2a, 2c
Table S1: List of diamonds used in this work, indicating the implantation energy, (column 2). Column 3: maximum implantation depth used in the modelling, , where it is assumed that the defects (N and NV) are uniformly distributed over the range , see justification in Sec. IV. In the main text, each sample is referred to in terms of the mean implantation depth, . Column 4: figures in which each sample appears.

i.2 Surface treatments

To form the H-terminated channels measured in main text Figs. 2b,c and 3, the diamond sample was first subject to a soft hydrogen plasma (7 minutes at 400 W, 10 Torr) optimised to make the diamond surface conductive via hydrogenation while preserving the integrity of the NV centres Hauf2011 (). The channels were then protected by a photoresist mask and the sample subject to a soft oxygen plasma (5 minutes at 50 W, 14 MHz RF with a 0.7 Torr O/Ar pressure) optimised to render the surface non-conductive via oxygen termination, while minimising the amount of etching of the diamond. No topography step was observed in atomic force microscopy (AFM) indicating an etching by the oxygen plasma below our AFM resolution (i.e.  nm). The resulting surface terminations are referred to as hydrogen-terminated (H) and oxygen-terminated (O) in the main text.

In main text Fig. 2d, the effect of other variants of oxygen termination are investigated. Namely, starting with a diamond (sample #2) oxygen terminated via an oxygen plasma as described previously (labelled as ‘Plasma’ in main text Fig. 2d), we applied the following steps: (i) acid cleaning (15 minutes in a boiling mixture of sulphuric acid and sodium nitrate); (ii) annealing at C in oxygen-rich atmosphere similar to the process used in Ref. Lovchinsky2016 (); (iii) cleaning in a piranha solution (mixture of 4 ml of sulphuric acid and 2 ml of hydrogen peroxide heated to 90C) for 10 minutes. The resulting surfaces are labelled in main text Fig. 2d as ‘Acid’, ‘O burn’ and ’Piranha’, respectively.

i.3 Device fabrication

The devices such as the one imaged in main text Fig. 3 were fabricated as follows. A stack of Ti/Pt/Au (thickness 10/10/70 nm) was evaporated onto the diamond (diamond #2, see Table S1) masked by a photoresist pattern. After lift-off of the photoresist leaving Ti/Pt/Au contacts (Fig. S4a), the sample was annealed at 600C for 20 mins in hydrogen gas (10 Torr). At such temperature, Ti atoms are able to diffuse into the diamond, and conversely carbon atoms are able to diffuse into the Ti layer, thus forming a TiC layer extending about 15 nm into the diamond (Fig. S4b, and Fig. S13 for the depth measurement). The thick Au layer serves as the primary contact material for electrical interfacing, while the Pt layer is introduced to act as a barrier preventing diffusion of Au and Ti atoms across the Ti/Au interface. The interest of such a process is to form a clean quasi-one-dimensional interface with the H-terminated channels subsequently fabricated (Fig. S4c). Together with the high expected work function of TiC ( eV Kim2015 ()), this results in more consistent formation of Ohmic contacts than with conventional 2D interfaces of H-terminated diamond with Au Hauf2014 () or Ti Akhgar2016 (), for instance. After the TiC/Pt/Au contacts were formed, the H-terminated channels were made via the process outlined above, i.e. hydrogen plasma to H-terminate the bare diamond surface, photo-lithography to protect the channels, and oxygen plasma to O-terminate the unprotected diamond surface. Finally, large Cr/Au (10/70 nm) contact pads partially overlapped with the TiC/Pt/Au contacts were evaporated onto the device for physical wire bonding.

Figure S4: Schematics illustrating the fabrication process: a Ti/Pt/Au stack is evaporated on the diamond (a) and annealed to form a TiC layer extending into the diamond (b), providing a low-resistivity interface with the conductive 2DHG formed by H-terminating the diamond surface with a hydrogen plasma (c).

i.4 Imaging set-up

The NV imaging set-up is a custom-built wide-field fluorescence microscope similar to that used in Refs. Simpson2016 (); Tetienne2017 (). The diamonds were glued on a glass cover slip patterned with a microwave waveguide (Fig. S5b), connected to a printed circuit board (PCB, see Fig. S5a) with silver epoxy. The diamond devices were electrically connected to the cover slip via wire bonding, and to the PCB board with silver epoxy. The voltage through the device under study (in main text Fig. 3) was applied using a source-meter unit (Keithley SMU 2450) operated in constant voltage mode. All measurements were performed in an ambient environment at room temperature, under a bias magnetic field generated using a permanent magnet (visible in Fig. S5a).

Optical excitation from a 532 nm Verdi laser was gated using an acousto-optic modulator (AA Opto-Electronic MQ180-A0,25-VIS), beam expanded (5x) and focused using a wide-field lens ( mm) to the back aperture of an oil immersion objective lens (Nikon CFI S Fluor 40x, NA = 1.3). The photoluminescence (PL) from the NV centres is separated from the excitation light with a dichroic mirror and filtered using a bandpass filter before being imaged using a tube lens ( mm) onto a sCMOS camera (Andor Zyla 5.5-W USB3). Microwave excitation was provided by a signal generator (Rohde & Schwarz SMBV100A) gated using the built-in IQ modulation and amplified (Mini-Circuits ZHL-16W-43+) before being sent to the PCB (see Fig. S5a). A pulse pattern generator (SpinCore PulseBlasterESR-PRO 500 MHz) was used to gate the excitation laser and microwaves and to synchronise the image acquisition. The total CW laser power at the sample was 300 mW, which corresponds to a maximum power density of about 5 kW/cm given the beam diameter. The optically detected magnetic resonance (ODMR) spectra of the NV layer were obtained by sweeping the microwave frequency while repeating the following sequence: s laser pulse, s wait time,  ns microwave pulse (except in Fig. S12 where it was varied); with total acquisition times of several hours typically.

Figure S5: (a) Photograph showing the PCB with microwave and DC inputs, as well as the permanent magnet used to apply the bias magnetic field . (b) Photograph of the diamond as mounted. The red box indicates a typical area used for imaging.

Ii Electric field mapping with NV ensembles

In this section, we describe a method to measure the full vector components of the magnetic and electric fields simultaneously, using ensembles of NV centres. To this aim, we first discuss the Hamiltonian of a single NV centre, then show how one can make use of multiple NV centres in a single crystal for maximal sensitivity to electric fields, and finally we describe our fitting and analysis procedure.

ii.1 Electrometry with a single NV

The spin Hamiltonian of the NV electron spin () in the presence of a magnetic field and electric field is given by


where are the spin-1 operators,  MHz is the temperature-dependent zero-field splitting,  GHz/T is the isotropic gyromagnetic ratio, and Hz cmV and Hz cmV are the electric susceptibility parameters Doherty2012 (); Doherty2013 (). Here is the reference frame of the NV defect structure as defined in Ref. Doherty2014 (), where is the major symmetry axis defined by the direction joining the nitrogen and the vacancy (which is along a crystal direction), and is a minor symmetry axis defined as being orthogonal to and also contained within one of the three reflection planes.

The frequencies of the and spin transitions, , can be computed by numerically diagonalising the above Hamiltonian. To facilitate the discussion, in what follows we will make use of the approximation


where the terms are defined as


This approximation is valid under the situation where Doherty2014 (), which is a good approximation for the cases explored in this paper.

As apparent from the expression of , sensitivity to the transverse electric field () is maximised when the longitudinal component of the magnetic field () is minimised. To a lesser extent, sensitivity to is also maximised when the transverse magnetic field () is minimised, due to the term in Eq. S3. To illustrate this, a simulation of the shift caused by an electric field of strength  kV/cm is plotted as a function of in Fig. S6a. Here the electric field is taken to be along the lab-frame axis as defined in Fig. S5b, i.e. the [001] direction of the diamond crystal, to mimic the electric field associated with surface band bending. Fig. S6a shows that there is indeed a dramatic decrease in the expected NV frequency shift with , e.g. from  MHz at zero magnetic field to  MHz with just  G, and a milder effect of the transverse magnetic field, for instance the shift is still  MHz under  G (with ). There is also a dependence on the orientation of in the transverse () plane as shown in Fig. S6b, which is related to the direction of both and relative to the defect minor axis , as captured by the term in Eq. S3. For transverse magnetic field strengths of the order of 50 G, the loss of sensitivity due to a non-optimal angle is relatively mild (a factor 2 at most), and as such careful alignment of in the transverse plane to match a given direction of is not critical. We stress that the shift induced by the longitudinal electric field () is usually much smaller than that from because , even though it doesn’t decrease with the application of a magnetic field.

Figure S6: Shift of the NV spin transition, , caused by an electric field of strength  kV/cm, plotted as a function of with in (a), and as a function of with in (b-d), where is the reference frame associated with the NV defect structure, with being along a direction and the minor axis Doherty2014 (). is indicated above each graph in terms of its Cartesian components in the lab frame as defined in Fig. S5b, i.e. is aligned along in (a,b), in (c) and in (d).

ii.2 Electrometry with NV ensembles

Clearly, measuring the two spin transition frequencies from a single NV centre is not sufficient to infer the six vector components of and simultaneously, without prior knowledge or additional measurements. On the other hand, using an ensemble of NV centres in a single crystal provides four independent measurements because can be along one of the four crystal directions, denoted as . As a result, with eight different frequencies measured there is enough information in principle to infer the six unknown field components, in addition to the zero-field splitting (). Experimentally, this requires all eight frequencies to be spectrally resolvable in the ODMR spectrum, which is achieved by applying a carefully aligned external magnetic field .

Precisely, is chosen to satisfy several criteria: (i) the direction of is chosen perpendicular to one of the directions (e.g., ) in order to maximise the sensitivity of the corresponding NV centres (family NV) to electric fields; (ii) the direction of within the said transverse plane is varied so that the projections of along the four NV axes are as distinct to each other as possible; (iii) the amplitude is chosen as a trade-off between electric field sensitivity which prescribes to be minimised (see Fig. S6) and sufficient spacing between adjacent ODMR lines (so they can be resolved given their linewidth). To illustrate the last two points, we calculated the minimum splitting between any of the 8 ODMR lines, with frequencies , defined by


as a function of when (i.e., ) and , as shown in Fig. S7a. In our experiments, we aimed for a minimum separation  MHz nominally, which allows variations in the ODMR frequencies across the field of view of several MHz to be measured (due to variations in the electric or magnetic field generated by the sample or the external magnet). As can be seen in Fig. S7a, this requires a magnetic field strength of the order of  G, with a wide range of directions allowed within the plane. A typical experimental ODMR spectrum taken in such conditions is shown in Fig. S7b where  G, the Cartesian components being expressed in the lab frame defined in Fig. S5b. Fig. S7b also defines our convention to label the 8 NV resonance frequencies and the 4 NV families. With this , the symmetry axes for the 4 NV families have unit vectors expressed as follows in the lab frame:

Figure S7: (a) Minimum separation between the 8 ODMR lines as calculated from Eq. S5, as a function of the transverse magnetic field relative to , which is the symmetry axis of the NV family optimised for electric field sensitivity (i.e., ). The red circle indicates the regime used in the experiments, which corresponds to a magnetic field with Cartesian coordinates  G in the lab frame . (b) Typical ODMR spectrum measured under this magnetic field , and conventions for the NV transition frequencies and for the NV families . The solid line is a fit with 8 Lorentzian lines with free amplitudes and widths.
Figure S8: (a) PL image of sample #2 showing a H-terminated channel between two TiC/Pt/Au contacts. (b,c) ODMR spectra averaged over the regions delimited by the square boxes in (a), with matching colours, recorded under a bias magnetic field  G, a voltage  V, and with the laser illumination spot positioned as in main text Fig. 3f. (d) ODMR spectra of the same regions as in (c) but under a magnetic field  G. In (b-d), the solid lines are data fits with 8 Lorentzian lines with free amplitudes and widths.

To illustrate how the presence of electric fields affects the ODMR spectrum, we show ODMR data corresponding to the largest electric fields observed in this work, allowing one to readily visualise the induced frequency shifts. We thus consider the situation examined in main text Fig. 3f, where a voltage  V was applied to the H-terminated device resulting in an electric-field feature near the drain contact. The PL image of the device is shown in Fig. S8a, and ODMR spectra of selected regions are shown in Figs. S8b-d under  G (b,c), which is the field optimised for simultaneous vector electrometry and magnetometry and used throughout the paper, and under  G for comparison (d). The two spectra in Fig. S8b are taken outside the H-terminated channel and show no visible shift in the NV frequencies (relative to each other), illustrating the good uniformity of over these length scales (m separate the two plotted regions). Under the conductive channel, however, there is a clear difference between the two spectra (Fig. S8c), where the blue spectrum corresponds to the large electric field seen in main text Fig. 3f. The central lines ( and , family NV) exhibit shifts by up to 4 MHz (for , based on a Lorentzian fit of the whole spectrum) as well as an extra broadening, compared with the reference red spectrum. The other lines are also slightly shifted, which is clearly observable for and (family NV), but more subtle for the other four lines (NV and NV). Under a non-optimised (Fig. S8d) where all the NV families have a significant component, the shifts are smaller, although they are still visible for NV. These observations are consistent with a change in electric field between the two regions, which shifts the NV frequencies to an extent that depends on the angle formed between and the NV symmetry axis (see Fig. S6a).

Figure S9: (a) ODMR spectrum obtained by computing the 32 transition frequencies (see details in text) and applying a Lorentzian lineshape with a fixed arbitrary amplitude and a full width at half maximum (FWHM) of 5 MHz (broad lines, comparable to those in the experiment) or 50 KHz (narrow lines, for ease of visualisation). The magnetic field is taken to be  G as in the experiment, while the electric field is null for the red spectrum, and  kV/cm for the blue spectrum. (b-d) Close-up views of (a).

The broadening of and seen in Fig. S8c is attributed to a combination of two effects. First, the electric field associated with surface band bending is expected to be non-uniform across the thickness of the NV layer (see Fig. S17b), causing an inhomogeneous broadening. The second source of apparent broadening is caused by a splitting of each NV line into two lines corresponding to the two sub-groups of NV centres distinguished by their orientation, i.e. N-V (where the vacancy is closer to the diamond surface than the nitrogen atom) or V-N (the vacancy is closer). In the absence of electric field, the NV transition frequencies are invariant under this inversion, but this is no longer the case in the presence of an electric field Doherty2014 (); Michl2014 (), especially when the axial component of the magnetic field is vanishing (). To illustrate this, we computed the transition frequencies for the 8 possible NV orientations (4 directions + NV/VN inversion) and generated an ODMR spectrum by applying a Lorentzian lineshape to each resonance. For completeness we also included the hyperfine interaction of the NV electron spin with the N nuclear spin (spin-), so that the total Hamiltonian of the system is

where is the nuclear spin operator,  MHz/T is the nuclear gyromagnetic ratio, and MHz and MHz are the hyperfine parameters Felton2009 (). One thus has a total of 32 lines in the ODMR spectrum (2 electron spin transitions for each of the 8 possible NV orientations and the 2 possible nuclear spin states), although some of them are nearly degenerate resulting in only 8 lines being usually resolvable under our experimental conditions (due to the intrinsic linewidth of 1-2  MHz in our samples Tetienne2018 () and additional power broadening). Illustrative simulated spectra are shown in Fig. S9a, obtained using  G as in the experiment, and an electric field either null (red spectrum) or of  kV/cm (blue). As expected, the presence of an electric field affects especially the central lines (called and according to our previous definition), which split further apart from each other, and additionally split into two sub-lines corresponding to the two possible orientations within family NV (i.e. N-V vs. V-N, see Ref. Doherty2014 (); Michl2014 ()), separated by nearly 9 MHz. A higher resolution spectrum (sharp lines in Fig. S9) reveal additional splittings by  kHz caused by the hyperfine interaction, which is highly suppressed for a purely transverse magnetic field. The other lines are also shifted overall by the electric field, where the high-resolution spectrum (Figs. S9c,d) reveals a small splitting induced by the orientation inversion (N-V vs. V-N) on top of the usual hyperfine splitting of  MHz.

Looking at the experimental spectrum (blue curve in Fig. S8c), we remark the presence of a shoulder on line which may be the signature of this inversion-asymmetry-induced splitting. However, the presence of strong electric field gradients mentioned above acts as a source of broadening which prevents detailed comparison with the theory. Furthermore, while the simulated spectra assumed a constant amplitude for each resonance, this is not the case in reality. Indeed, the transition strength (Rabi frequency) depends on the initial and final states of the NV spin system and on the direction of the microwave field relative to the NV defect structure, and can therefore vary significantly between different transitions and in the presence of an electric field, especially for NV centres with a small axial magnetic field (e.g. NV). For instance, this explains why the two lines and exhibit different amplitudes even in the low electric field case (see red spectrum in Fig. S8c, we verified that had a smaller Rabi frequency than , by a factor , under the same microwave power). This complicates the analysis of the experimental spectrum when in the presence of an a priori unknown electric field distribution. Consequently, for imaging purpose we simply fit the experimental spectrum with a sum of 8 Lorentzian lines (solid lines in Figs. S8b-d) with free frequencies , amplitudes and widths, thus providing a mean value for each electron spin transition incorporating the effect of hyperfine and inversion-asymmetry splittings. The 8 frequencies are then used to deduce the magnetic and electric field vector components, as will be detailed in the next section.

ii.3 Fitting procedure

Figure S10: (a) Maps of the frequencies obtained by fitting the ODMR spectrum at each pixel. The region shown is the same as that in Fig. S8a, and the bias magnetic field is  G. (b,c) Maps of the magnetic field components (b) and out-of-plane electric field component (c) obtained after reconstruction. In (b) a plane subtraction was applied to remove the background field .

To infer the unknown values from the measured frequencies , we seek to minimise the root-mean-square error function


where are the calculated frequencies obtained by averaging over both the NV orientation for the corresponding NV family (i.e., the axis is ), and over the nuclear spin state (), that is,


Each frequency is obtained by projecting and into the NV reference frame knowing the orientation of the axis (), numerically computing the eigenvalues of the Hamiltonian (II.2), and deducing the electron spin transition frequencies.

For the fitting (i.e. minimisation of ), we used a bound constrained optimisation based on the fminsearch function in Matlab. While it is possible in principle to fit all 7 parameters, we found that the fit would not reliably converge for all the pixels within the field of view, resulting in very noisy electric field images (hence large uncertainty). We conclude that the noise level of our measurements (the uncertainty in the data is 20-60 kHz for a single pixel, see Table S2) is insufficient to perform a full vector fit for each pixel. On the other hand, we found that fixing and letting only free, which is the component expected from surface band bending, solved the convergence issues resulting in much smoother images. In all the data shown in the main text, we therefore fit only the five parameters , while fixing . We also note that because the measurement and model are averaging over the N-V and V-N sub-families, there is an ambiguity in the overall sign of the determined electric field. In other words, and yield the same ODMR spectrum, hence one measures when . Solutions to determine the sign of include using preferentially-oriented NV centres Michl2014 (); Lesik2014 () and applying an external electric field Dolde2011 (). As an illustration of the reconstruction process, Fig. S10a shows the 8 frequency maps obtained for the same area as in Fig. S8a. In this example, a voltage  V was applied to the device, corresponding to a current A flowing through the conductive channel. This current produces a non-trivial magnetic field distribution through the Bio-Savart law, resulting in a complex landscape in the frequency maps. After reconstruction following the procedure outlined above, one obtains the maps of the 3 vector components of (Fig. S10b) as well as (Fig. S10c). The magnetic field maps are consistent with the magnetic field expected from a uniform current density in a flat wire Tetienne2017 (). On the other hand, the map shows a very different distribution with minimal cross-talk with the magnetic field, illustrating the effectiveness of the reconstruction method.

ii.4 Error analysis

(kHz) (kHz) (kHz)
#1 630 45 23
#2 380 60 22
#3 570 63 47
#4 270 90 53
#5 220 61 20
Table S2: Residual error from fitting the measured ODMR frequencies without including the Stark effect (second column) and when including a normal-to-the-surface electric field (third column). The fourth column indicates the standard error on the from fitting the ODMR spectrum with Lorentzian functions. The data shown corresponds to the measurements described in main text Fig. 2a, i.e. O-terminated diamond implanted at different depths, and all values are averaged over a large number of pixels.

We now discuss the sources of error in our reconstruction method. First, it is interesting to examine the residue after fitting. This is summarised in Table S2 for the data shown in main text Fig. 2a, which gives the case where the Stark effect is not included in the Hamiltonian, i.e. , and the case used in the paper where the Stark effect is included but the electric field is constrained to be perpendicular to the diamond surface, i.e. . Also indicated in the uncertainty (standard error) in determining the frequencies , denoted as , obtained from the Lorentzian fit to the ODMR spectrum and averaged over the eight lines. Clearly, the data is not well fit without including the Stark effect, with being 5-30 times larger than the measurement uncertainty . Instead, when one includes the Stark effect ( only), one reduces the residue by up to an order of magnitude, which is clear evidence that measurable electric fields are present in our sample. The residue is still a factor 2-3 larger than , which can be due to one or a combination of the following effects. First, the electric field might be not exactly along the axis, e.g. due to surface roughness, diamond miscut, or non-uniform density of surface or bulk defects, although fitting and as well as did not seem to reduce the residue . Second, other corrections to the Hamiltonian may need to be considered, e.g. the effect of mechanical strain in the diamond (see discussion in Sec. II.5). Third, there may be systematic errors in the measured , for instance due to each ODMR line comprising several sub-lines with uneven amplitudes (they are assumed even in our analysis, see Fig. S9 and corresponding discussion) or due to partial overlapping between adjacent lines (leading to systematic errors if the actual line shape is not exactly Lorentzian as it is assumed, which can arise e.g. from the electric field gradient experienced by the NV ensemble). Nevertheless, the fact that simply adding a single parameter () to the fit function reduces the residue by an order of magnitude suggests that this parameter captures very well the essence of the problem.

We listed above the possible causes for systematic errors in the measurements. On top of that, there is a statistical error in the measured that originates from photon count noise in the ODMR spectrum, which is usually close to the photon shot noise limit in our experiments but can also be affected by readout noise and dark counts of the camera depending on the exact conditions. To characterise this noise, we simply calculate the standard deviation from a large ensemble of pixels, which is found to be in the range 20-60 kHz in our experiments (see Table S2). This uncertainty on the translates into an uncertainty on . Again, we characterise this by taking the standard deviation of the map. Doing so on a small region gives an uncertainty of the order of 1 kV/cm, i.e. less than 1 % of the mean value, corresponding to the statistical (pixel-to-pixel) noise. By evaluating the standard deviation over larger areas, one captures statistical noise as well as actual spatial variations of , due e.g. to spatial variations in the density of surface defects. Such variations can be seen in main text Fig. 2b and are of the order of 10 kV/cm, i.e. a few percents of the mean value. This is the uncertainty quoted in the main text and shown as error bars in the main text figures.

Finally, we discuss the consequences of the Stark effect from built-in electric fields on the accuracy of magnetometry measurements. In our measurements, because is treated as a fit parameter, neglecting the Stark effect leads to a biased estimation of in order to minimise . Taking sample #1 from Table S2 as an example, we obtain when including , against when neglecting the Stark effect. This shows that neglecting the Stark effect may lead to systematic errors of the order of T in magnetometry measurements, motivating the need for precise characterisation of built-in electric fields for high-precision magnetometry applications, and minimising the sensitivity to electric fields by careful alignment of an external bias magnetic field.

ii.5 On the effect of strain

Mechanical strain acts as an effective electric field in the NV spin Hamiltonian Doherty2012 (), which simply adds to the true electric field and therefore may in principle contribute to the electric fields measured in this work. In this section, we estimate the effect of strain in our samples induced by bulk defects, surface defects as well as dilational dislocations, and conclude that it is negligible compared to the electric field due to surface band bending.

ii.5.1 Strain induced by bulk point defects

In the simplest picture, a point defect is an isotropic point source of expansion or contraction of the lattice Stoneham1968 (). In continuum mechanics, the displacement field of this expansion/contraction is analogous to the electric field of a point charge and takes the form Stoneham1968 ()


where is known as the defect strength and is the position of the defect. This displacement field is the solution of the elastic equation


Since the displacement field induced by a point defect is irrotational , the displacement field can be expressed in terms of a scalar potential, , via . In which case, the elasticity equation becomes (see references Love2011 () for background)


Note the parallel between and the electric potential.

The defect strength is related to the change in volume induced by the defect Stoneham1968 ()


where is the Poisson ratio of the solid. The change in volume can be obtained from X-ray crystallography that determines the solid’s lattice parameter as a function of the defect concentration or by ab initio calculations.

The components of the strain field are defined with respect to the displacement field via


where denotes the -vector component of the displacement field and denotes the -coordinate. The stress is calculated from the strain via


where and are Lame’s constant and the shear modulus, respectively.

Assuming that the displacements due to different defects are sufficiently small that they may be added in superposition, the elasticity equation for a continuous number density of defects is simply


Consider now a semi-infinite slab in the plane, extending vertically from to . If the slab has a uniform density of defects, the solutions of the elasticity equation immediately follow as


All other strain and stress components are zero. Thus, the defects produce uniform stress and strain that is principally in the vertical direction (assuming ) and proportional to the density of defects.

For the N defect in diamond, using the x-ray data ppm Lang1991 () and the diamond Poisson ratio . Using  GPa and  GPa, the stress components as a function of concentration (units: GPa/ppm) are and repectively. Since the spin-stress susceptibility parameters are on the order of 1 MHz/GPa, and the N concentrations in our implanted samples are on the order of 10 ppm, we get that the NV frequency shifts induced are  kHz, equivalent to  kV/cm in terms of effective electric field. This is smaller than our measurement uncertainty and therefore negligible. The strengths of other common defects (N, NV and NV) are smaller than N, as predicted by the ab initio calculations in Biktagirov2017 (). Thus, even considering multiple defect species, it is unlikely that strain will induce a significant shift in the spin resonances.

ii.5.2 Strain induced by surface defects

To model surface defects, reconsider the slab but with a defect density that decays exponentially from the surface


In this case, the elasticity solutions are


Again, the principal stress is in the vertical direction and is proportional to the density of defects. As the density decays exponentially from the surface, defects densities that are highly localised to the surface do not generate stress inside the solid, and their effect on the NV spin resonances can therefore be neglected.

ii.5.3 Stress induced by dilational dislocations

A dislocation that induces pure compression/dilation is equivalent to a line of point defects with uniform density. So, if we say that the dislocations have a uniform density in the transverse plane of a slab, this is equivalent to a uniform plane of point defects. Since we have concluded that the stress/strain at a given height in the slab depends only on the local density of defects, and we expect dislocations induced by polishing/surface damage to be localised to the surface, then we can conclude that dislocations do not influence the NVs sufficiently below the surface. Thus, regardless of their strength, dislocations may also be ignored.

Iii Supplementary data

iii.1 Electrical characterisation of the devices

Figure S11: Electrical characterisation of a representative H-terminated channel between two TiC/Pt/Au contacts. (a) Capacitance under a bias voltage  V, as a function of the frequency of the AC signal. (b) Capacitance at 20 Hz as a function of the bias voltage (laser on). (c) Source-drain current as a function of the applied DC voltage.

The electrical properties of the hydrogen terminated channels were characterised first by DC measurements to determine the device resistance (Keithley 2400 SMU), and then with a small AC voltage test signal ( mV peak to peak) applied by a LCR meter (HP 4284A) to determine the complex device impedance (Z). The resistance of the devices was modified by the  nm laser illumination required by NV measurements. Upon illumination the resistance of all devices dropped from  k to  k ( V). Green laser illumination is known to optically pump charge carriers in nitrogen doped diamond through deep defects and surface defects Nesladek1998 (). The photo-induced carriers in our devices were observed to have a long lifetime with photo-currents taking minutes to decay. Prolonged laser exposure resulted in a gradual and permanent increase in the device resistance. After multiple days of continuous laser illumination the devices typically exhibited dark resistances of  M (illuminated resistances of  k), which can be explained by a decrease in dark hole density with no change in the photo-induced carrier density. The reduction of hole density indicates a loss of surface acceptors either by laser induced removal of the hydrogen termination or by the formation of a surface contamination that prohibits transfer doping of the diamond.

The TiC contacts to the hydrogen terminated ribbon exhibited an Ohmic response within the applied DC bias range of to  V shown in Fig. S11c. The sub linear response seen above  V is consistent with a drift saturation velocity of  cm s Hauf2014 (). Despite these stable Ohmic contacts, AC measurements revealed a non-zero parallel capacitance in the diamond device. At low frequencies the device exhibits a large negative capacitance (Fig. S11a, red line), which becomes positive ( pF) above  kHz. We interpret this capacitance as a result of the geometrical differences between the TiC contacts, which extend  nm into the diamond, and the conductive H-terminated channel, which is confined within  nm Hauf2014 () of the surface, resulting in charge build up at subsurface TiC-diamond interface, below the 2D hole gas. The frequency response of the negative capacitance is inconsistent with a parasitic inductance, and is understood to result from the transient relaxation current () after a small voltage step not decreasing monotonically (i.e. at some time, , after the voltage step). The longer times required before results in the measured negative capacitance disappearing with high frequency measurements, as the voltage changes are quick enough to ensure  Ershov1998 (); Jonscher1986 (). The microscopic cause of this unusual transient response is device specific, for our devices we postulate that the negative capacitance is a result of carrier trapping dynamics across the metal TiC and p-type diamond interface similar to the cause of negative capacitance seen in silicon based p-n junctions Lemmi1999 ().

Fig. S11b shows the low frequency negative capacitance response to the applied DC bias. The symmetry around suggests that TiC-diamond interfaces at both contacts behave consistently under a particular bias direction. Electric field images of the biased interfaces (main text Fig. 3) suggest that this negative capacitance corresponds with the increased electric field and occurs when the interface is forward biased, i.e. when the p-type H-terminated channel is positively biased with respect to the TiC. This indicates that either holes are injected into the TiC, electrons are injected into the diamond or both. As TiC is typically metallic Leroy2006 (), these holes would rapidly recombine, however electrons injected into the conduction band far away from the diamond surface may have a longer lifetime. Trapping of these electrons in ionised donor defects would re-distribute charge within the diamond resulting in a change in the electric field at the NV layer and cause a positive transient response in the current due to the hole density increasing to maintain charge neutrality with the adsorbed acceptor layer. Laser illumination would pump any filled trap states reducing the measured electric field, and increase the magnitude of the positive component of the current transient which results in the increase in the negative capacitance at low frequencies shown in Fig. S11a (blue line).

iii.2 Electric field vs laser intensity

To assess the effect of the laser used for the NV measurements on the band bending (through photo-ionisation of bulk defects), we measured the electric field as a function of laser intensity for diamond #2 (with no applied voltage). Precisely, we kept the peak laser power (300 mW) and laser pulse duration (s) constant and varied the dark time , which varied the laser duty cycle and hence the average power (Fig. S12a) according to


Since photo-currents have been measured to be long lived (minutes compared to the microsecond time scale of the laser pulsing sequence), we expect varying the average laser power via pulsing to be equivalent to varying the laser power in a CW experiment. We observed a decrease in by about 15% when the duty cycle was increased from 2.5% to nearly 100%, i.e. a factor 40 increase in average laser intensity. This indicates that the laser has a measurable effect on the band bending. All measurements reported in the main text were performed with a duty cycle close to 1 in order to maximise the signal-to-noise ratio, and as such they include a small but measurable effect of the laser. For simplicity, photo-ionisation effects are ignored in our modelling (section IV), which can be a source of discrepancy when comparing to the experiment.

Figure S12: The effect of the laser on the measured electric field. (a) Pulse sequence for ODMR defining the laser pulse time and the dark time. (b) Electric field measurements for oxygen (blue) and hydrogen (orange) terminated diamond for different laser duty cycles, solid line is an exponential fit to the data.

iii.3 Electric field vs diamond etching

Here we investigate the effect of etching the diamond on the electric field. Namely, we used a diamond implanted at  nm (sample #2) and applied two steps of etching ( nm per step) over partly overlapping patches.

The first etching step was done by removing the TiC/Pt/Au contacts off sample #2 via acid etching (15 minutes in a boiling mixture of sulphuric acid and sodium nitrate to remove the Pt/Au layers, followed by minutes ultrasonication in :: NHOH : HO : HO to remove the TiC), as illustrated in Fig. S13a. Following the etching, the diamond was optically clear indicating that there was no metal left, and a recess of  nm was measured by AFM in the regions formerly occupied by the TiC/Pt/Au stack (Figs. S13b,c).

Figure S13: (a) Illustration of the etching process: the TiC layer formed by annealing is etched away, leaving a recess. (b) Atomic force microscopy (AFM) image of a region thus etched. (c) Line cut extracted from (b) revealing a  nm etching step.

The second etching step was done by reactive ion etching (RIE) through a Au mask, which was subsequently removed via acid etching (15 minutes in a boiling mixture of sulphuric acid and sodium nitrate). An etching step of  nm was measured by AFM and optical profilometry. The etching mask was chosen to partly overlap with the regions etched during the first etching step, so that the diamond eventually had regions that were etched only once ( nm, due to either step 1 or step 2) or etched twice ( nm, due to both steps) as well as non-etched regions (Fig. S14a). Because an acid cleaning was applied as the final step of the process, all the surfaces should be oxygen terminated, yet with possibly different density of surface defects.

Figure S14: (a) Topography image of the diamond surface obtained by optical profilometry, highlighting (A) a non-etched region, (B) regions etched once and (C) a region etched twice. (b) Electric field map of the same area as in (a). The plotted range is capped to 280 kV/cm to highlight spatial variations, with actual values reaching up to 320 kV/cm. (c) Corresponding PL image. (d) Line cuts extracted from (a-c), taken along the dashed lines shown in the images.

The topography, electric field, and PL images of a selected area are shown in Figs. S14a-c, respectively, with line cuts shown in Fig. S14d. The electric field is found to be correlated with the surface height, namely the more the diamond was etched, the smaller the electric field became. On the other hand, the PL shows a less trivial correlation, where it is larger in the regions etched once ( nm etching) than in the regions etched twice (30 nm) or not etched at all. To interpret these results, let us first consider the effect of the first etching step, which decreased the electric field by and increased the PL by a factor , despite etching away 15 nm off the estimated depth range of the NV centres,  nm. This implies a dramatic reduction in band bending, such that for instance the depletion region for NV decreased from 27 nm (i.e. NV only at depths  nm, an 8-nm-thick layer) to just 5 nm (NV at  nm, where  nm has been etched). The measured is averaged over the NV depth distribution, which is deeper before the etch hence probes a flatter part of the band bending, leading to only a small decrease in (by ) upon etching even though is expected to be significantly smaller (by a factor of with our example numbers). A possible explanation for such a dramatic reduction in band bending invokes the presence of sub-surface defects present in the non-etched region and induced by the hydrogen plasma initially applied to the sample, which may have resulted in diffusion of hydrogen atoms into the diamond Stacey2012 (). The etching would then remove these defects to leave only the surface-termination-induced defects as the source of band bending.

Assuming the second etch results in the same surface as the first etch, with the same density of surface defects, we would expect the PL to decrease relative to the first etch, which it does by a factor . The assumptions made above (5 nm depletion layer and NV depth range  nm before etching) are clearly over simplistic as they predict no PL at all after a 30 nm etch. Most likely, the implantation profile (at energy 10 keV for this sample) produced a tail that extends deeper than 35 nm. The further decrease in upon the second etch can be explained by an increase in the band bending depth (i.e., the length scale of the bending) while the total bending (i.e. the energy offset at the surface) remains unchanged. This is expected as at this stage a significant number of nitrogen defects (donors) has been etched away, leaving a diamond with a lower nitrogen concentration hence a more dilute space charge density. While further work is needed to fully understand the multiple competing effects, these results illustrate the value of in-situ electric field measurements, which provide new insight into band bending.

Iv Modelling of band bending at the diamond surface

In this section, we first present our model of the oxygen- and hydrogen-terminated diamond surfaces, then describe how band bending and electric field are calculated within this model, and finally discuss some results.

iv.1 Models of the diamond surface

Figure S15: (a) Band diagram of a nitrogen implanted diamond near the surface, showing the conduction band minimum (CBM), valence band maximum (VBM), the N charging level, the NV charging level, the NV charging level, and the constant Fermi level . To the left of the vertical axis are: a depiction of the density of states associated with the sp surface defects, which have a total areal density , and an ionised density ; a depiction of the adsorbed acceptor molecules present close to the VBM for hydrogen-terminated diamond, with an ionised areal density (), and inducing a 2D hole gas in the diamond. (b) A cartoon depicting the surface and nitrogen induced defects considered in the band modelling (see text).

iv.1.1 Oxygen-terminated diamond

The surface of materials often have electronic energy levels that are distinctly different from the bulk material bands. In semiconductors, these surface states can exist within the band gap of the material, as depicted in Fig. S15. Any unoccupied (occupied) surface states (with areal density ) below (above) the bulk Fermi level will ionise. However, to maintain charge neutrality within the crystal the total charge of the ionised surface states needs to be compensated by an opposite charge from the bulk. In order for charge to build up in the bulk the bands bend towards the surface, which simultaneously reduces the density of the ionised surface states, until charge neutrality is reached. Typically the Fermi level at the surface is pinned close to the surface state energy unless a density of donors that fully compensate the density of surface states can be introduced into the material. In the case of oxygen-terminated diamond we observed a strong electric field which we ascribe to an upward band bending indicative of a high density of surface states. This is commensurate with published electrical measurements of Fermi level pinning at diamond interfaces Zheng2001 (); Garino2009 (); Mori1991 ().

Information about the density and energy of any surface states present on oxygen terminated surfaces will then be required to model the expected electric field. The cross section of a (001) diamond surface consists of carbon bonded in a zig-zag pattern between the first and second atomic layer. Each surface carbon atom has two dangling bonds that are satisfied by bonding to an oxygen termination layer. Oxygen can terminate the diamond in a number of ways: an oxygen atom can bond to two surface carbon atoms in an ether like arrangement (depicted in Fig. S15b), it can double bond to a single surface carbon atom in a ketone like arrangement or a hydroxyl group can terminate a single surface carbon atom. Each of these different arrangements introduce different surface states into the diamond band gap, yet only ether-like terminations are expected to result in upward band bending in nitrogen doped diamond (where the Fermi level is close to the nitrogen donor level at VBM + 3.75 eV) with an unoccupied state at VBM + 3.4 eV Sque2006 (). Such a high unoccupied energy level is unable to explain the large band bending that was measured on all oxygen terminated samples and thus we don’t consider the effect of any oxygen induced surface states on the diamond electronic structure. Recently a primal defect, the sp defect, of the diamond surface has been identified to have a much lower acceptor energy around VBM + 1.5 eV Stacey2018 (). The fundamental structure of this defect is shown in Fig. S15b, and occurs when a surface carbon atom is removed from the lattice. The dangling bonds that remain then form double bonds with the exposed sub-surface carbon. As the energy of this defect state is well below the nitrogen donor level, and the density is expected to be large,  cm, this defect makes a good candidate to explain the measured electric field under oxygen terminated surfaces. The actual energy of the sp defect surface state is dependent on its local environment, and DFT modelling predict energies within a range from  eV to  eV. To capture this range we consider the defect density of states (DOS) to be Gaussian around  eV with a full width at half maximum of  eV, in agreement with the shape of the resonance peak seen in X-ray absorption experiments probing this defect Stacey2018 ().

iv.1.2 Bulk defects

The surface charge is compensated by the dielectric response of the diamond, , the ionisation of any defects within the band gap, predominantly nitrogen donors, and any charge carriers (which is practically zero until the Fermi level is within 1 eV of a diamond band). The donors within the diamond are formed by the implantation of nitrogen ions, which has a fixed areal density (fluence) of nitrogen cm in these experiments. For shallow implants with an implant energy below  keV, the depth dependent nitrogen density profile is in general non trivial Lehtinen2016 () and is particularly sensitive to the angle of incidence, which here is nominally from the surface normal but is partly randomised by the roughness of the diamond surface. We chose to approximate the nitrogen density profile by a rectangular function bounded by the surface on one end () and a maximum implant depth parameter, , on the other end. Such a simple function is a reasonably good approximation to both theoretically and experimentally determined profiles in our situation Toyli2010 (); Lehtinen2016 (); FavaroDeOliveira2015 (), and avoids the introduction of multiple parameters. The only parameter is , which can be related to the implant energy via . To maximise the match between the data presented in the main text and theory, a value of was chosen, hence a mean depth , also providing a reasonable agreement with previous studies Tetienne2018 (); Lehtinen2016 (); Toyli2010 (). The volumetric density of nitrogen in the diamond is then calculated as . The nitrogen implant forms a range of different defects within the diamond: substitutional nitrogen (N), nitrogen-vacancy centres (NV) and vacancy defects (which can be single, double or multi-vacancy chains). The majority of the implanted nitrogen ends up as N with only % turning into NV for high nitrogen fluences Pezzagna2010 (), resulting in densities = and = for the N and NV defects, respectively. The N acts as a deep acceptor state,  eV ( eV) away from the conduction band. Whilst the NV defect can become both positively and negatively charged, allowing to act as both a super deep donor ( eV) and mid gap acceptor ( eV) Deak2014 (). For simplicity we do not consider other vacancy related defects. Although the expected density of vacancy defects is larger than the density of nitrogen related defects after implantation, multiple annealing steps at C and C were undertaken to minimise the quantity of vacancy defects and the final density is expected to be low Goss2005 (); Tetienne2018 (); DeOliveira2017 ().

iv.1.3 Hydrogen-terminated diamond

Having presented our model of the oxygen-terminated surface, we now move on to describe the hydrogen-terminated surface. The clean and ordered hydrogen-terminated diamond surface forms a strong dipole field that reduces the diamond’s work function, moving the vacuum energy level below the diamond conduction band (observed as a negative electron affinity Takeuchi2005 ()). An observed result of this low work function is that it is energetically favourable for electrons in the diamond valence band to be excited to adsorbed acceptor molecules in atmospheric conditions. The total negative charge of ionised adsorbed acceptors (), depicted in Fig. S15a, needs to be compensated by a positive charge within the diamond, and much like under an oxygen-terminated surface the diamond bands bend upwards towards the surface. As this bending is not pinned to an energy level in the diamond band gap, it often extends all the way into the diamond valence band, forming a two-dimensional hole gas confined within  nm of the surface RgenRistein2006 (); Hauf2014 () and resulting in the observed surface conductivity of hydrogen-terminated diamonds Koslowski2001 (). In intrinsic diamond, holes in the valence band are the only positive bulk charge that provides a significant contribution to compensation of the surface charge. In nitrogen doped diamond, however, the ionised nitrogen defects will also partially compensate any surface charge. With a high enough density of nitrogen, the Fermi level at the surface will move far enough from the valence band that no hole layer is formed and the diamond surface ceases to be conductive Grotz2012a (). In our samples, however, the primal sp defects should still be present because the hydrogenation is unlikely to remove them Stacey2018 (). As a consequence, the sp defects will once again pin the surface Fermi level position. Like with oxygen termination, the total charge in the acceptor layer is a function of the Fermi level at the diamond surface Newell2016 (), however to simplify dealing with the effect of both surface acceptors and ionised sp defects, we consider to be constant. We can relate the measured hydrogen ribbon resistivity () to the density of holes at the diamond surface, , using , where  m and  m are the hydrogen ribbon length and width for our devices. The conductivity of the device is given by where is the hole mobility in the ribbon. Hole mobilites between  cmVs are commonly reported Sato2013c (); Pakes2014 (). Device resistances of our ribbon ranged from initial values of  k to  M after extensive measurement, not taking into account any contact resistance. This provides an estimate for the 2D hole density in the range of  cm. Accounting for the  keV implant in sample #2 ( nm), numerical integration returns that is within the range to  cm.

In all cases the upward bending bands will affect the response of the implanted NV defects, de-ionizing the near surface NV where the band bending is strong Schreyvogel2016 (), and creating a space charge layer of ionised N (the green circles in Fig. S15). The space charge layer creates a significant electric field around the NV defects, it is this field that can be observed as ODMR shifts. To estimate the amount of NV and the size of the electric field that the NV experience, the shape of this bending needs to be calculated.

iv.2 Calculating band bending

Poisson’s equation provides a classical definition of how a charge distribution changes the electric potential within a dielectric material:


where is the vacuum permittivity and is the relative permittivity of the material. Solving this equation determines how the bands bend given a certain surface potential (or equivalently, a normalised potential , related to via Eq. S39), i.e. the boundary condition at the diamond surface. The shape of the bending is only determined by the charges within the diamond, i.e. if one solves for and redefine the surface at discarding the higher potentials, the solution is the same as solving for  lüth2014solid (). We take advantage of this and over-solve all our bands with a large then redefine the surface at the appropriate potential to compensate a given areal density of charge at the diamond surface.

To determine the charge distribution in the diamond we need to consider the densities of electrons, holes and any ionised defect states (negatively charged acceptors and positively charged donors). The density of electrons, holes and ionised defects is determined by the filling of electronic states in the diamond valence and conduction bands and the defect states, respectively. The filling of electronic states is calculated by convolving the energy dependent DOS with the probability of a state being occupied by an electron at a given energy. The probability of a state being occupied (unoccupied) is given by the Fermi-Dirac distribution function,


and its complement,


where is the Boltzmann constant, is the temperature, and the functions are centred around the Fermi level () of a material. The Fermi level can be determined by requiring that the material be charge neutral at the Fermi level, giving


where () is the density of holes (electrons), is the density of ionised donors, and is the density of ionised acceptors in the bulk of diamond in thermal equilibrium. We can calculate the electron and hole densities using the effective mass () approximation for the valence and conduction bands DOS StreetmanBanerjee ():