# Influence of (N,H)-terminated surfaces on stability, hyperfine structure,
and zero-field splitting of NV centers in diamond

Wolfgang Körner<EMAIL_ADDRESS>Fraunhofer Institute for
Mechanics of Materials IWM, Wöhlerstraße 11, 79108 Freiburg, Germany  Reyhaneh
Ghassemizadeh  Fraunhofer Institute for Mechanics of Materials IWM,
Wöhlerstraße 11, 79108 Freiburg, Germany  Daniel F. Urban  Fraunhofer
Institute for Mechanics of Materials IWM, Wöhlerstraße 11, 79108 Freiburg,
Germany  Christian Elsässer  Fraunhofer Institute for Mechanics of Materials
IWM, Wöhlerstraße 11, 79108 Freiburg, Germany  University of Freiburg,
Freiburg Materials Research Center (FMF), Stefan-Meier-Straße 21, 79104
Freiburg, Germany

###### Abstract

We present a density functional theory analysis of the negatively charged
nitrogen-vacancy (NV-) defect complex in diamond located in the vicinity of
(111)- or (100)-oriented surfaces with mixed (N,H)-terminations. We assess the
stability and electronic properties of the NV- center and study their
dependence on the H:N ratio of the surface termination. The formation energy,
the electronic density of states, the hyperfine structure and zero-field
splitting parameters of an NV- center are analyzed as function of its distance
and orientation to the surface. We find stable NV- centers with bulk-like
properties at distances of at least $\sim 8$ Å from the surface provided that
the surface termination consists of at least 25% substitutional nitrogen
atoms. The studied surface terminations have a minor effect on the ground
state properties whereas the NV orientation has major effects. Our results
indicate that axial NV centers near a flat 100% N-terminated (111) surface are
the optimal choice for NV-based quantum sensing applications as they are the
least influenced by the proximity of the surface.

###### pacs:

67.30.er, 07.55.Ge, 71.70.-d

##  I Introduction

The negatively charged nitrogen-vacancy center NV- in diamond is a point-
defect complex with excellent potential for the use in spatial-atomic-
resolution quantum magnetometry [1, 2, 3] and in solid-state-based quantum
computing.[4] It consists of a substitutional nitrogen atom in the diamond
crystal structure with a vacant nearest-neighbor carbon site and an additional
electron. The focus of this study is on NV- centers as functional elements in
scanning magnetic-field sensors. A detailed review on the various sensing
mechanisms using NV-centers can be found in Ref. [5] In principle, for sensing
desirable to bring them as close as possible to the diamond surface in order
to be able to position this atomic magnetic-field probe as close as possible
to the external magnetic field to be measured. Unfortunately, NV- centers may
loose their negative charge state near surfaces with negative electron
affinity (EA). A prominent example is the case of a hydrogen-terminated
diamond surface where the EA is between -1 and -1.3 eV according to
experiments.[7, 8, 6, 9, 10] Besides the negative EA, additional surface
states may lie deeply in the band gap and interact with the electronic levels
of the NV- center.

Various theoretical studies[11, 12, 13] have determined surface states for H,
OH, F, and N terminations for (001), (111) and (113) surfaces. In summary, all
the studied hydrogenated surfaces show deep levels below the conduction band
that disturb the levels of the NV- center, whereas fluorinated and
nitrogenated surfaces have only few surface states and highly positive EA. For
a comprehensive overview on the surface states see Figure 1 of Ref. [11]. A
summary of the EA can be found in Table 1 of Ref. [14].

Our study deals with surfaces terminated by H atoms, N atoms, or mixtures of
these two elements. Nitrogen atoms are most promising as they create the least
amount of surface states while the detrimental hydrogen atoms can hardly be
avoided in any experimental diamond-growth process. By considering surface
terminations with different N:H ratios one can estimate to what extent H may
be tolerable when aiming at the experimental realization of NV- centers very
close to a surface. Therefore we want to answer the following questions: (i)
How much hydrogen can be tolerated in a mixed (N,H)-termination of the surface
in order to keep the negative charge state of the NV- center? (ii) At what
distance from the surface does a NV- center have to be situated so that its
functional properties are only marginally disturbed by the proximity of the
surface?

Our density functional theory (DFT) analysis is restricted to NV centers which
are at most about 15 Å below the surface due to computational costs. Thus we
investigate the case of very shallow NV centers which has recently already
been realized experimentally. Müller et al.[15] reported single NV centers at
20-30 Å and Ofori-Okai et al.[16] reported even implantation depths of $<$10
Å. In the following the NV centers are investigated at all possible sites
relative to the (001) and (111) surfaces in atomistic supercell models. By
evaluating their total energies it is clarified whether certain sites near a
surface are preferred with respect to sites in the bulk crystal or not.
Furthermore, the electronic density of states (DOS), the hyperfine structure
(HFS) and the zero-field splitting (ZFS) parameters are evaluated. With these
results we quantify the influence of the surfaces, and we assess whether the
electronic level structure of the NV- center near the surface is still intact.

The manuscript is organized as follows: In Sec. II the details of the DFT
calculations and the atomistic supercell models are described. The results of
the NV- stabilization including EA, total energy and DOS calculations are
presented in Secs. III.1, III.2, and III.3, respectively. The results for the
ZFS calculations are reported and discussed in Sec. III.4, and the HFS
analysis is presented in Sec. III.5. Section IV summarizes our findings. The
Appendix gives a brief compilation of important properties of the tensors of
hyperfine structure parameters $\boldsymbol{A}_{ij}^{I}$ and zero-field
splitting parameters $\boldsymbol{D}_{ij}$.

##  II Theoretical approach

###  II.1 Supercell models

In this study the two most prominent cases of diamond surfaces are considered,
namely the (001) and (111) surface orientations. As in our previous work[17]
we build the atomistic supercell models from the bulk primitive cubic cell of
diamond with a lattice constant of $a$= 3.567 Å which agrees with the
experimentally determined one of Holloway et al..[18] Cubic diamond (Ramsdell
notation[19] 3C) has a layered structure with an ABC stacking sequence of
carbon double layers in a [111] direction. A (111) oriented hexagonal bulk
unit cell contains 6 carbon atoms. The supercell models for the (111)
orientation contain $6\times 6$ hexagonal unit cells in the a-b plane with
lattice vectors of length 15.13 Å. The supercell dimension in [111]-direction
(c axis) is 40.2 Å which includes an additional vacuum range of 10–12 Å,
depending on the considered surface termination. For 100% N-termination we
consider 13 double layers of carbon atoms and an additional monolayer of
nitrogen on each of the two surfaces in the supercell, cf. Fig. 1(a). This
model then comprises 1008 atoms in total, 936 C and 72 N atoms. For 100%
H-coverage the model contains 13 double layers plus an additional terminating
carbon layer with H atoms attached in positions “on top” of C atoms, cf. Fig.
1(e). This model then contains a total of 1080 atoms, thereof 72 H atoms. For
mixed (N,H) coverage, the terminating layer consists of a mixture of nitrogen
atoms and C-H pairs, as illustrated in Fig. 1(c). The size of the slab model
is about 26 Å for 100% N termination and 28 Å for 100% H termination. We
enumerate the inner 13 carbon double layers as indicated by the black numbers
on the left of Fig. 1(a). For clarity, only half of the supercell model is
displayed.

(a)(b)(c)(d)(e)(f) Figure 1: (Color online) Atomic structure models of the (a)
(111)-oriented surface and (b) (001)-oriented surface with 100% N termination.
The terminating atomic structures for a 50% mixed ratio of H and N as well as
for the 100% H termination are illustrated in (c)–(f). Yellow spheres
represent carbon atoms, small red spheres hydrogen atoms, and blue spheres
nitrogen atoms. The position of a carbon vacancy is indicated by a large gray
sphere. In order to indicate the positions of the NV centers relative to the
surface, the single layers for (001) and the double layers for (111) are
enumerated. For the (001) surface the seven different (symmetry inequivalent)
positions of NV centers are indicated, which differ in their relative
positions with respect to the N atoms on the surface. For the (111) surface we
plot the example of an axially oriented NV center in the fifth double layer,
denoted NV5(ax), as well as a NV center with opposite orientation in the
fourth double layer, denoted by VN4(ax). NV2(ba) denotes a NV in basal
orientation in the second double layer.

For the (001) orientation we use supercells with $6\times 6$ cubic primitive
units in the a-b plane. The corresponding lattice vectors have a length of
15.13 Å. The supercell dimension in c direction amounts to 44.3 Å including a
vacuum range of about 10–12 Å. The diamond supercell with 100% N termination
was constructed with 33 carbon monolayers plus an additional layer of nitrogen
on each of the two sides of the slab. This gives 72 N and 1188 C atoms. Fig.
1(b) displays half of the supercell with the enumeration of the single C
layers on the right side, in order to indicate and specify the various
positions of NV centers. For 100% hydrogen coverage the 33 central carbon
monolayers are terminated by an additional carbon layer with H atoms attached
on top of C atoms, cf. Fig. 1 (f). This supercell contains 1260 C and 72 H
atoms in the supercell summing up to 1332 atoms in total. In the case of mixed
(N,H) we follow the surface construction of Refs.-[39, 40] suggesting the
formation of N-C-H bonds at the surface, cf. Fig. 1(d). Moreover, we examine
homogeneous and clustered distributions of N atoms on the surface, which is
discussed in detail in Sec.VI.1. The lowest energy configurations of the mixed
surfaces, being those with a homogeneous distribution of N-C-H groups, are
illustrated in the inset of Fig. 3.

Note that for the pure hydrogen or nitrogen terminations, our atomic surfaces
structures correspond to those already discussed in the literature (eg. see
Ref. [12] for (111) and Ref. [20] for (001)).

(a)(b) Figure 2: (Color online) (a) Top view on a specific (low energy)
structure model of a (111) surface with a coverage of 50% H and 50% N. An
axial and a basal orientation of a NV center are shown as a black circle and
ellipse, respectively. (b) Top view on the atomic structure model of the (001)
surface with 100% nitrogen termination. The seven different possible locations
of NV relative to the surface are illustrated. The possibilities shown here
repeat every four monolayers. However, the differences between the variants
become less important with growing distance from the surface.

In order to determine the interplay of the surface and the NV- center we
considered all symmetry-inequivalent NV positions within the supercells. These
positions vary in distance and orientation with respect to the surface. For
the (111) orientation and mono-atomic termination there are four possible
orientations due to the high symmetry of the surface. NV centers are either
oriented parallel to the [111] direction (called axial) or inclined by about
109.5∘ with respect to the c-axis (called basal). In both axial and basal
configurations there are two possibilities to place the NV center, with either
the substitutional N atom or the C vacancy closer to the surface. These are
denoted by NV and VN, respectively. The axial NV have the full C3v symmetry as
in the bulk crystal while the basal NV has a reduced C1h symmetry due to the
presence of the surface. In the case of mixed (H,N) terminations all the NV
centers have the reduced C1h symmetry which implies a splitting of certain
energy levels near the surfaces.

For the (001)-oriented surface with 100% nitrogen termination the situation is
more complicated. There are seven sets of NV centers which differ in their
relative orientations with respect to the N-N pairs on the surface, as
illustrated in Fig. 2(b). Due to the layer stacking sequence in c direction,
there are equivalent orientations every four layers. For surfaces with mixed
(N,H)-termination the situation becomes even more complicated. Of course, with
increasing distance to the surface all the properties of the various NV
centers assimilate to those of a NV center in the bulk crystal.

In the following sections we will analyze the properties of NV- centers in
terms of their distance to the surface. For this study we define the distance
between a NV center and the surface by the distance of the vacancy
(representing the geometrical center of the defect complex) to the first full
layer of carbon atoms. For example, the NV center denoted VN3I in Fig. 1a) has
a distance of 1.78 Å and NV2(ba) in Fig. 1b) has a distance of 4.12 Å. In our
supercell models the NV centers can be placed at a maximum distance of about
14 Å to the surface. The distance between two periodic images of NV centers in
the a-b plane is about $15.1$ Å, i.e. the length of the respective lattice
vectors of the supercells.

###  II.2 Computational details

For the structural relaxation of the atomistic supercell models and the
calculation of the physical parameters of interest we use the Vienna Ab Initio
Simulation Package (VASP)[22, 21]. The Bloch waves of the valence electrons
are expanded in a plane-waves basis (with a cutoff energy of 420 eV) and the
interactions of the valence electrons and the ionic cores are included by
projector-augmented-waves (PAW) potentials.[23] The exchange-correlation
energy and potential are treated in the generalized gradient approximation
(GGA), as given by Perdew, Burke and Ernzerhof (PBE).[24]

For all the supercells the Brillouin-zone integrals are evaluated only at the
$\Gamma$-point with a Gaussian broadening of the energy levels by 0.05 eV. The
positions of the atoms in the constant supercell volume were relaxed until the
residual forces acting on them were less than 0.03 eV/Å and the energy
difference between two consecutive ionic relaxation steps was less than
$10^{-5}$ eV.

The computation of the HFS tensor components $\boldsymbol{A}_{ij}^{I}$ and the
ZFS tensor components $\boldsymbol{D}_{ij}$ was done using subroutines
implemented in VASP. A summary of the key formulae can be found in the
Appendix.

For the calculation of the electron affinity $\chi$ we used the pure surface-
slab models without NV centers. To first order, $\chi$ is given by

$\chi=E_{vac}-E_{CBM},$ (1)

where $E_{vac}$ is the vacuum level and $E_{CBM}$ is the energy of the
conduction-band minimum (CBM). $E_{vac}$ is determined by calculating the
average electrostatic potential perpendicular to the surface and by
subsequently taking the plateau value in the vacuum region (see e.g. Fig. S4
of the supplemental material of Ref. [11].) $E_{CBM}$ is obtained from the
local density of states (DOS) of carbon atoms in the respective bulk regions
of our models. Since the PBE functional underestimates the band gap of diamond
by more than 1 eV, a common practice is to add the experimental band gap value
of 5.47 eV[25] to the calculated valence band maximum.[26, 27] Theoretical
results using this correction scheme are indicated by the label _+corr_ in
Table 1. Our results do not include any such post-processing correction. A
detailed discussion on the calculation of the EA can be found in Ref. [26].

In order to obtain and model a charged NV- center one can either consider
charged supercells that contain an additional electron or include a second
substitutional nitrogen atom in the structure. The advantage of the latter
approach is that the slab in total is charge neutral. However, the extra N
atom introduces additional symmetry breaking which will affect the surfaces
and more importantly the NV-. In our extensive supercell models of about 1000
atoms we could not find a site that is distant enough such that the N atom is
not affecting the properties of the NV-. We have tested both approaches and
obtained qualitatively very similar results for the formation energies, DOS
and hyperfine constants but the transversal component E of the ZFS differs
substantially since it is very sensitive to differences in the charge
densities of the three next-neighbor C atoms to the vacancy of the NV-. We
have chosen the charged supercell approach for this study and can therefore at
least describe the experimentally important case of axial NV centers at (111)
surfaces adequately (see Sec. III.4.2). The drawback of this methodology is
that we have to deal with a charged slab which manifests itself in the
presence of small slopes in the ZFS parameters D and E as well as in the HFS
constants Aii, as function of the NV center position. However, the respective
changes are smaller than 1% across the inner slab, so that the quantities
appear constant in the bulk regions (see Figs. in the following Secs. III.4
and III.5). A detailed discussion of several approaches to simulate charged
diamond-surface slabs can be found in Ref. [28].

##  III Results and discussion

###  III.1 Electron affinity

Termination | (001) H |  (111) H |  (001) N | (111) N
---|---|---|---|---
This work, PBE | -1.02 | -1.0 | 4.0 | 3.5
Experiments | -1.30a | -1.27b |  |
LSDA+corrc | -1.96 | -2.01 |  |
HSE06d | -1.7 | -1.6 | 3.5 | 3.2
Table 1: Comparison of electron affinities $\chi$ for (001) and (111) surfaces
with either H- or N-termination. The experimental data is taken from aRef. [7]
and bRef. [25], the theoretical literature data stems from cRef. [27] and
dRef. [28].

Figure 3: (Color online) Electron affinity for (001) and (111) surfaces with
mixed (N,H) coverage as function of nitrogen (hydrogen) proportion. The three
insets exemplary show the considered (001) surface terminations. The dotted
line at 0.3 eV marks the average offset of our results compared to
experimental values obtained for pure H coverage, cf. Tab. 1.

We evaluated the EA for mixed (N,H)-coverage as function of the partial
nitrogen concentration (or the N:H ratio) of the surface termination using the
supercell models without any added point-defect complex. Pure hydrogen
coverage leads to a negative EA and pure nitrogen coverage to a positive EA
(see results in Table 1). Our point of interest is the N:H ratio at which the
EA changes sign. Therefore, we have investigated NV- centers in slab models
with 16.7%, 20%, 25%, 50% and 100% nitrogen on the surface. Figure 3 displays
our results for the two considered surface orientations (001) and (111). The
sign change occurs at approx. 33% nitrogen for the (111) surface and at
approx. 25% for the (001) surface.

This critical value for the N:H ratio was also observed in our subsequent
simulation of the NV-. For slab models with less than 25% nitrogen on the
surface the NV- defect complex was not stable even at the maximum distance of
about 14 Å from the surface. The NV center could not bind the extra electron
and, instead, the negative charge accumulates in the vacuum region of the
supercell. However, this is not a real physical effect but an artefact of the
periodic supercells which do not allow for a complete withdrawal of the
charge. For 25% nitrogen or more we obtained self-consistent solutions with
the extra electron bound to the NV defect complex.

Chou et al.[28, 12] calculated the EA using the HSE hybrid functional and
reported the sign change of the EA at 45% nitrogen and 55% hydrogen on the
(111) surface. Their results for mono-atomic coverage are included in Tab. 1
for comparison. Our EA results for pure H coverage differ from the
experimental values only by an offset of approximately 0.3 eV. Shifting our
results (that were obtained using the PBE functional) by 0.3 eV yields the
sign change of the EA at approx. 40% nitrogen for the (111) surface and at
approx. 33% for the (001) surface, as illustrated by the dotted black line in
Fig. 3.

Experimental evidence for shallow stable NV- centers for slightly negative EA
were recently reported.[29] The experimental results for the EA are in
approximate quantitative agreement with the results of our simulations which
include the NV center and an extra electron in the supercell. The discussion
in the following sections is restricted to details of stable NV- centers near
surfaces with N:H ratios of 1:3, 1:1 and 1:0, corresponding to 25%, 50% and
100% nitrogen, respectively.

###  III.2 Formation energies of NV- centers

Figure 4: (Color online) Formation energies of NV defect complexes with
different orientations in the vicinity of (N,H)-terminated surfaces relative
to their formation energies inside diamond at 14 Å distance to the surfaces.
Top: (111)-oriented surface with data points starting from double layer 1;
Bottom: (001)-oriented surface with data points starting from single layer 3.
For positions closer to the surfaces no stable NV center-like defect complexes
were obtained.

Figure 5: (Color online) Formation energies of NV centers in the vicinity of
100% N-terminated (001) surface. The data is divided into sets according to
the different orientation with respect to the N atoms on the surface, cf. Fig.
2(b). The gray dashed-dotted line connects the lowest energies of the sets as
function of distance (which yields the red line in the lower panel of Fig. 4).

We investigated the influence of the proximity of the surface on the energetic
stability of the NV center by a comparison of total energies. Figure 4
displays the formation energies of NV- centers as function of their distance
and orientation with respect to the surface. For every distance considered,
only the value of the lowest-energy orientation is plotted. Values are given
relative to the formation energy of NV- at the innermost bulklike layers.
Here, the double layer $n=7$ is the reference for the (111) surface, while the
monolayer $n=16$ is the reference for the (001) surface (cf. Fig. 1).

NV- centers that are incorporated closer than approx. 8 Å to the surface
partially vary significantly in their formation energies by more than $\Delta
E=\pm 0.1$ eV. Close proximity to the (111) surface causes a strong surface
reconstruction which may yield large changes in energy. For example, NV1(ax)
has $\Delta E$ $\approx$ 2.5 eV. The typical geometric arrangement of the NV-
defect complex is not conserved in this case, and thus it does not make sense
to classify it as a NV- center. Correspondingly, for the (001) surface, we
found that defects in layers 1 and 2 lead to substantial rearrangement of the
surrounding atoms, so that they can no longer be classified as a NV- defect.
Therefore, we only show data points for which the local geometry of a NV-
complex is preserved.

For the (001) surfaces with NV-orientation we display the data for 25%, 50%
and 100% nitrogen concentrations at the surface in the lower panel of Fig. 4.
It can be seen that the formation energy is practically independent of the
surface chemistry. The same behavior is observed for the VN-orientation as
well as for the axial and basal NVs for the (111) surfaces. As all these
results are very similar, we refrain from plotting them for all considered N:H
ratios for the sake of clarity.

The alternating up and down in energies observed between 2 and 6 Å from the
surface in the case of (001)-orientation can be understood from the seven
different sets of data, corresponding to the different possible orientations
of the NV- with respect to the N atom pairs on the surface, as illustrated in
Fig. 5. Due to the specific surface reconstruction in (001)-orientation NV-
centers can be accommodated in different relaxation modes within the same
monolayer or with respect to the monolayers above or below. For instance, at a
distance of 4.5 Å the NV5I is much lower in energy than NV5II and the NV3I and
NV5I are energetically more favorable than NV4 and NV6I in the neighboring
monolayers. However, as can be seen, the differences between the sets of data
become less pronounced with growing distance from the surface. This emphasizes
the short-range effect of the surface reconstruction which does not disturb
stable isolated NV- centers which are approx. 8 Å away from the surface.

We conclude that NV- centers can be stabilized very close to the (001) and
(111) surfaces since neither towards the surfaces nor towards the bulk region
there is a significant energy gradient. Especially, the positions NV3I and
NV5I for (001) and NV1(ba) for (111) are the most stable positions. However,
very close to the surfaces the electronic levels of the NV- centers are
strongly modified as will be shown in the following.

###  III.3 Eletronic levels

In bulk diamond, far from any surfaces or other defects, NV- centers have
characteristic sharp electronic levels lying deep within the band gap. The
single-electron levels $a_{1}(2)$, $e_{x}$ and $e_{y}$ are well separated from
the band edges.[30, 12] Near surfaces the situation changes since surface
states appear. For example in the case of hydrogen termination these lead to a
wide electronic band below the CBM of bulk diamond.[11] Such bands interfere
with $e_{x}$ and $e_{y}$ if the NV- center is spatially located close enough
to the surface.

Our (111) surface models with mixed termination and at least 25% nitrogen have
an open band gap like in bulk diamond and deep lying NV- defect states that
are neither close to the conduction band nor the valence band. This was
already reported for pure nitrogen termination by Chou et al.[12] for the
(111) surface. For our (001) slab models, surface states appear at approx. 0.8
eV below the CBM for mixed termination with 25% nitrogen. Their position
changes to 0.3 eV below the CBM for 50% nitrogen.

Despite the open band gap, a NV- center that is close ($<7$ Å) to a surface
causes structural relaxations that alter its electronic levels. In the bulk
interior, the $e_{x}$ and $e_{y}$ levels are degenerate since they are formed
by the equivalent three carbon atoms next to the vacancy.[30] This degeneracy
of $e_{x}$ and $e_{y}$ is lifted near the surfaces due to the reduced
symmetry. Figure 6 shows this energy splitting of the single-electron levels
$e_{x}$ and $e_{y}$ for NV- centers near the (001) surface with 100%
N-termination. A maximal splitting of 250 meV of the single-electron levels
$e_{x}$ and $e_{y}$ was obtained at position NV3 decreasing to 50 meV at
position NV9. For NV- centers at even larger distance to the surface, we
observe that all of them also show this minimal level splitting of $\approx
50$ meV as a consequence of using a finite slab model geometry for the
simulation. NV- centers with basal orientation near a (111) surface behave
similarly like NV- centers near the (001) surface. By contrast, the threefold
symmetry is retained for NV- centers with axial orientation near a perfectly
flat infinite (111) surface with pure N termination.

We refrain from giving absolute values for the electronic levels not only
because it is well known that the PBE functional is already underestimating
the band gap of bulk diamond by more than 1 eV but also because there are
further problems. In finite slab models of a few nm size there are confinement
effects, self interaction effects in lateral direction, and also the presence
of two polarized surfaces due to NV- in the model structure. A thorough
discussion of these effects for the NV center in diamond was presented by Chou
et al.[28] In Fig. 5 of their paper the authors also present the defect levels
of NV- in the vicinity of various surfaces which agree qualitatively with our
results.

(a)(b) Figure 6: (Color online) (a) Splitting of the single-electron levels
$e_{x}$ and $e_{y}$ for NV- centers near the (001) surface with 100% N-
termination. The single-electron levels $a_{1}(2)$ are only marginally
affected by the vicinity of the surface and therefore not shown. (b) Level
structure of the NV- with values taken from Ref. [30]. NV- energies are given
with respect to the VBM.

###  III.4 Zero field splitting

####  III.4.1 Longitudinal ZFS component D

The zero-field splitting component $D$ of the singlet and triplet states of
the ${}^{3}A_{2}$ ground state of a NV- center located in a bulk diamond
crystal was experimentally measured as $D$=$\frac{3}{2}D_{zz}$= 2.872(2) GHz
with high precision by Felton et. al.[31]. Our theoretical values for bulk
diamond exceed the measured value by about 3 to 6% depending on the supercell
size and other computational details. A discussion of the limitations and
approximations can be found in our previous work[17] and in Ref. [32] which
addresses the central problem of the so-called spin contamination.

In Fig. 7 the axial component $D$ of the ZFS for our surface models normalized
to the value $D_{0}$ obtained from the respective bulk supercells is shown.
Near the surface the $D$ values are reduced, where the reduction is in all
cases more pronounced for the VN-oriented NV- centers than for the NV-oriented
ones. This can be understood from the fact that the main contribution to the
ZFS originates from the three C atoms next to the vacancy. Each of them has a
magnetic moment of about 0.52 $\mu_{B}$, amounting together to 78% of 2
$\mu_{B}$. For the NV- centers with NV-orientation these three C atoms lie
further away from the surface and are thus situated in a more bulk-like
neighborhood than the ones of the VN-orientation.

The influence of the surface on the axial component $D$ of the ZFS is limited
to about 8 Å below the surface irrespective of the surface orientation and its
chemical termination. As shown in Fig. 7, a very weak dependence on the N:H
ratio is found for both (001) and (111) surfaces.

Figure 7: (Color online) Longitudinal component $D$ of the ZFS of a NV- center
for different positions relative to the surface. Top panel: (111)-orientation;
Bottom panel: (001)-orientation.

####  III.4.2 Transversal ZFS component E

In perfect diamond and zero external magnetic field $B_{{\rm ext}}$ along the
NV- axis the $m_{{\rm S}}=\pm 1$ states are degenerate, where $m_{S}$ denotes
the spin projection on the NV- axis. In experiments a magnetic field can be
sensed directly by measuring the Zeeman-splitting (namely by a splitting of
5.6 MHz per 1 Gauss field) of the $m_{{\rm S}}=\pm 1$ states.[5] However,
crystal imperfections may lead to finite values of the transversal component
$E$ and thus additional splittings as discussed e.g. in our previous work.[17]
The total splitting of the energy levels for the spin $S=1$ system reads

$\displaystyle E_{m_{{\rm S}}=\pm 1}$ $\displaystyle=$ $\displaystyle
D\pm\sqrt{E^{2}+(g\mu_{B}B_{{\rm ext}})^{2}}$ (2)

and $E_{{\rm m_{S}=0}}=0$.

Figure 8: (Color online) Transversal component E of the ZFS of a NV- center
with axial orientation near (111)-oriented surface.

The calculation of $E$ is computationally demanding. It does not only suffer
from systematic errors like the aforementioned spin contamination[32] but in
our situation the charged NV- defect induces dipole charges at the two
surfaces of the slab models. Theses dipoles cause an artificial electric field
across the slab. Only for axial NV- centers the effect of the resulting dipole
field cancels due to symmetry, since all three C atoms near the vacancy, which
are responsible for the main contribution to $E$, experience the same field.

Therefore, we present the calculated values of the transversal component $E$
of the ZFS for the ${}^{3}A_{2}$ ground state only for NV- centers with axial
orientation. In Fig. 8 one can see that $E$ depends crucially on the specific
(111) termination. For 100% N-termination the threefold symmetry is conserved
and thus $E=0$ MHz, independent of the distance to the surface. For mixed
(N,H)-terminations this symmetry is broken which leads to small splittings on
the order of 10 MHz close to the surface. However, these values quickly decay
to 0 MHz at larger distances to the (111) surface.

###  III.5 Hyperfine splitting connected to ${}^{14}{\rm N}$ and ${}^{13}{\rm
C}$

Calculated hyperfine structure constants for the NV- centers in the cubic
diamond crystal are in good accordance with recent experimental work of Felton
et al.[31], as discussed in our previous work[17]. This is a good starting
point for the study of the influence of the (111) and (001) surfaces on the
HFS. Among the NV-center properties investigated in this work, the HFS
constants are the least dependent on the proximity to the surfaces. As an
example, the HFS related to 13C for the practically most relevant case of
axial NV centers near (111) surfaces is shown Fig. 9. For distances of approx.
4 Å or more the NV- centers take their typical bulk values. For NV- centers
closer to the surface we observe that the Aii values are more affected for the
VN-oriented ones. This is conceivable since for VN-orientation the C atoms are
closer to the surface than in the NV-orientation. Basally oriented NV centers
show modified HFS for distances of approx. 6 Å, so one double layer deeper
into the bulk.

The results for the Hyperfine splitting connected to ${}^{14}{\rm N}$ are very
similar to the ones of ${}^{13}{\rm C}$, which means that the change of the
$A_{ii}$ is restricted to the first 6 Å. However, since the absolute
frequencies ($A_{11}$= -2.16 MHz and $A_{33}$= -1.73 MHz[17]) are two orders
of magnitude smaller compared to the $A_{ii}$ of 13C the relative increase is
in few case of order 20% very close to the surfaces ($A_{33}$= -1.35 MHz for
VN2(ax) near 100% nitrogen (111) surface).

We obtain a similar result for the (001) surface. Here the deviation of Aii is
5% for ${}^{13}{\rm C}$ and 10% for ${}^{14}{\rm N}$ for the NV- centers
closer than 5 Å to the surface, regardless of the surface termination. In
conclusion, we observe that the HFS constants converge quickly to their bulk
value with increasing distances to the surface.

Figure 9: (Color online) Hyperfine structure constants $A_{ii}$ for 13C for
axially oriented NV- centers near the (111) surfaces with different
termination.

###  III.6 Limitations of our model

The results presented above where obtained within the framework of DFT (PBE-
GGA) calculations and do not take into account some important dynamical
effects which are encountered in real diamond NV-center systems. In order to
achieve a high magnetic sensitivity, the coherence time (T2) of the NV center
needs to be long enough (in the range of $\mu s$). An important decoherence
mechanism which affects T2 is caused by fluctuating unpaired surface spins due
to impurities[33] and not atomically flat surfaces.[34] Romach et al.[35]
studied the surface induced noise and associated the high frequency part to
phonon effects because this is significantly suppressed at 10 K. They assigned
the temperature independent low frequency part of the noise to a surface spin
bath. Since in our static DFT calculations the atomic nuclei are kept at rest
($T=0$ K) and we restrict our investigations to atomically flat and defect-
free surfaces near the NV centers of interest, the dynamical effects mentioned
above are excluded. Nevertheless, our results can be considered as an ideally
static limit: a perfectly flat N-doped surface does not have unpaired spins
and thus avoids dynamical surface disturbances on NV centers.

##  IV Summary

We have studied the influence of the proximity of diamond (111) and (001)
surfaces with various chemical (N,H)-terminations on the formation energies,
densities of states, zero-field splittings and hyperfine structure constants
of NV- centers in diamond. For surfaces with less than approx. 25% nitrogen
the simulations did not yield stable negatively charged NV- centers but
neutral ones instead.

The analysis of formation energies shoes that NV- centers may be stable at
positions a few Å below the (001) and (111) surfaces. There is no steep energy
gradient, meaning the energy differences with respect to the bulk region are
very small. The result for the formation energies is that formation of NV-
centers closer than 6 Å to the surface occurs at a maximum energy gain of 0.15
eV. On the other hand, for distances larger than $\sim$ 8 Å the NV- formation
varies only in the range of 20 meV with respect to the bulk region.

Moreover, very close to the surfaces the electronic levels of NV- centers get
modified. The most prominent feature in the electronic structure of NV centers
is the splitting of the single-electron levels $e_{x}$ and $e_{y}$.

Furthermore, the HFS parameters $A_{ii}$ connected to ${}^{13}C$ (and
${}^{14}N$) result in a typical decrease (increase) of 5% to 20% relative to
the bulk values. These quantities are influenced by the proximity of the (001)
and (111) surfaces only in a very short range of approx. $4$ Å.

The ZFS parameters are affected by the surfaces on a longer range. The axial
component $D$, which is indicative for the singlet-triplet splitting, can be
reduced up to approx. 25% near the surfaces. For axially oriented NV centers
at (111) surfaces the transversal component $E$ of the ZFS is small (few MHz).

Finally, the dependence of all considered quantities on the surface
orientation (001) or (111) and its specific chemistry, namely the N:H ratio of
the surface termination, is weak.

In conclusion, we believe that our analysis improves the understanding of very
shallow NV centers in diamond and the data for the DOS, ZFS and the HFS may
help to interpret experimentally obtained spectra. Our results show that axial
NV centers near the flat 100% N-terminated (111) surface are the ideal choice
for NV-based quantum sensing applications since they feel the least influence
by the proximity of the surface due to their high symmetry.

##  V Acknowledgments

Financial support for this work was provided by the Fraunhofer Lighthouse
Project Quantum Magnetometry (QMag). The calculations were partially performed
on the supercomputer ForHLR funded by the Ministry of Science, Research and
the Arts Baden-Württemberg and by the Federal Ministry of Education and
Research.

##  VI Appendix

In this appendix we shortly discuss the influence of different surface
configurations at fixed N:H ratio and summarize the notation and calculation
of the HFS and the ZFS. For more detailed information on the HFS and ZFS in
diamond see e.g. Refs. [36, 37, 17].

###  VI.1 Properties of different (N,H)-configurations at the (001) surface

(a)(b)(c)(d)(e)(f)(g)(h) Figure 10: (Color online) (a)-(f) Different
configurations of N and H atoms on the (001) surface with a fraction of 25% N.
(g) Electron affinity (EA) and (h) energy difference relative to configuration
(a) which has the lowest-energy. Configuration (a) was chosen for the study
described in the main text.

In order to obtain the optimal surface configuration for cases of mixed
(N,H)-coverage we examined different orderings of N and H atoms on the (001)
surface with a fraction of 25% N atoms. As a starting point we considered
surfaces with homogeneously distributed N-C-H bonds and subsequently increased
the clustering of N atoms as depicted in Fig. 10(a)-(d). In addition, the
formation of N-N bonds as shown in Figs. 10(e) and (f) was tested. A general
trend in all these configurations is that with increasing clustering of N
atoms on the surface the formation energy increases. In the case of
configurations with N-N bonding the energy difference increases even up to 5
eV. Therefore, the N-C-H bonding on the (001) surface is energetically
preferred and was chosen for the further calculations in our study (see (a) of
Fig. 10).

Moreover, although the change of the EA for the tested surface configurations
is small, one can see that whenever the clustering of N atoms is increased
(cf. Figs.10(a)-(d)) also the electron affinity of the surface increases.
Note, that the arrangements (e) and (f) have direct N-N bonds which increase
the EA even more. In that sense, a purely N terminated surface can be regarded
as a huge cluster of N atoms and it is conceivable that such a configuration
has a positive electron affinity (see Fig. 3).

Throughout the paper we have used the energetically favorable highly ordered
(N,H)-configurations for the (001) and for the (111) surface slab cells with
NV defect. Note that already at a few carbon layers below the surface the
individual (N,H)-configuration becomes unimportant.

###  VI.2 Hyperfine interaction

The HFS tensor $\boldsymbol{A}^{I}$ describes the interaction between a
nuclear spin SI and the electronic spin distribution. The hyperfine
interaction between a nuclear spin SI (i.e. the nuclear spin of 13C or 14N)
and the electronic spin distribution Se (here NV- defect state) can be modeled
with the Hamiltonian $H_{{\rm HFS}}=S_{e}\boldsymbol{A}^{I}S_{I}$. The
hyperfine structure tensor components A${}_{ij}^{I}$ for a nucleus $I$ are

$A_{ij}^{I}=\frac{\mu_{0}\gamma_{I}\gamma_{e}}{2S}\int
d^{3}r~{}n_{S}({\bm{r}})\left[\left(\frac{8\pi\delta(r)}{3}\right)+\left(\frac{3x_{i}x_{j}}{r^{5}}-\frac{\delta_{ij}}{r^{3}}\right)\right],$
(3)

where the first term in brackets is the Fermi-contact term and the second
bracket is the magnetic dipole-dipole term. Here, $n_{S}$ denotes the spin
density associated with spin state S, $\mu_{0}$ is the vacuum magnetic
permeability, $\gamma_{e}$ the gyromagnetic ratio of the electron and
$\gamma_{I}$ the gyromagnetic ratio of the nucleus. In this work the values
$\gamma{(^{13}{\rm C})}/2\pi$ = 10.7084 Mhz/T and $\gamma{(^{14}{\rm
N})}/2\pi$ = 3.077 Mhz/T are used.[38]

###  VI.3 Zero-field-splitting

The ZFS describes the loss of degeneracy of the electronic levels of the spin
triplet state of the NV- (with different values of magnetic spin quantum
number mS= 0,+1,-1) in the absence of an external magnetic field. It arises
from the presence of unpaired electrons and can be modeled by the Hamiltonian

$H_{{\rm ZFS}}=\frac{\mu_{0}g^{2}\mu_{B}^{2}}{4\pi r^{5}}\bigl{[}3(S_{1}\cdot
r)(S_{2}\cdot r)-(S_{1}\cdot S_{2})r^{2}\bigr{]},$ (4)

which describes the effect of electron-electron repulsion by magnetic dipole-
dipole interaction. Here, $r=r_{1}-r_{2}$ is the spatial distance between the
spins, $S_{i}=\frac{1}{2}[\sigma_{x},\sigma_{y},\sigma_{z}]$ is the spin
operator vector of particle i, $\sigma_{j}$ (j=x, y, z) are the Pauli
matrices, and g is the Landé factor. One can separate the spatial and spin
dependencies in equation (4) and write the Hamiltonian in the form $H_{{\rm
ZFS}}=S\boldsymbol{D}S$, where $S=S_{1}+S_{2}$ is the total spin and
$\boldsymbol{D}$ describes the dipolar spin-spin interaction. The tensor
$\boldsymbol{D}$ is symmetric and traceless and thus can be diagonalized.

In general, one can split $H_{{\rm ZFS}}$ into a longitudinal component of the
magnetic dipole-dipole interaction (named $D$) and a transversal component
(named $E$):

$D=\frac{3}{2}D_{zz}~{}~{}{\rm and}~{}~{}E=\frac{1}{2}(D_{xx}-D_{yy}),$ (5)

where the $D_{ii}$ are the diagonal elements of the tensor $\boldsymbol{D}$.

## References

  * [1] G. Balasubramanian, I. Y. Chan, R. Kolesov, M. Al-Hmoud, J. Tisler, C. Shin, C. Kim, A. Wojcik, P. R. Hemmer, A. Krueger, T. Hanke, A. Leitenstorfer, R. Bratschitsch, F. Jelezko, and J. Wrachtrup, Nature 455, 648 (2008).
  * [2] J. R. Maze, P. L. Stanwix, J. S. Hodges, S. Hong, J. M. Taylor, P. Cappellaro, L. Jiang, M. V. G. Dutt, E. Togan, A. S. Zibrov, A. Yacoby, R. L. Walsworth, and M. D. Lukin, Nature 455, 644 (2008).
  * [3] V. M. Acosta, E. Bauch, M. P. Ledbetter, C. Santori, K.M. C. Fu, P. E. Barclay, R. G. Beausoleil, H. Linget, J. F. Roch, F. Treussart, S. Chemerisov, W. Gawlik, and D. Budker, Phys. Rev. B 80, 115202 (2009).
  * [4] V. Jacques, P. Neumann, J. Beck, M. Markham, D. Twitchen, J. Meijer, F. Kaiser, G. Balasubramanian, F. Jelezko and J. Wrachtrup, Phys. Rev. Lett. 102, 057403 (2009).
  * [5] M. Lee, J. Yoon and D. Lee, Magnetometers - Fundamentals and Applications of Magnetism, IntechOpen(2019) DOI: http://dx.doi.org/10.5772/intechopen.84204
  * [6] L. Rondin, G. Dantelle, A. Slablab, F. Grosshans, F. Treussart, P. Bergonzo, S. Perruchas, T. Gacoin, M. Chaigneau, H.-C. Chang, V. Jacques, and J.-F. Roch, Phys. Rev. B 82, 115449 (2010).
  * [7] F. Maier, J. Ristein, and L. Ley, Phys. Rev. B 64, 165411 (2001).
  * [8] K.-M. C. Fu, C. Santori, P. E. Barclay, and R. G. Beausoleil, Appl. Phys. Lett., 96, 121907 (2010).
  * [9] M. V. Hauf, B. Grotz, B. Naydenov, M. Dankerl, S. Pezzagna, J. Meijer, F. Jelezko, J. Wrachtrup, M. Stutzmann, F. Reinhard, and J. A. Garrido, Phys. Rev. B, 83, 081304(R) (2011).
  * [10] S. Kawai, H. Yamano, T. Sonoda, K. Kato, J. J. Buendia, T. Kageura, R. Fukuda, T. Okada, T. Tanii, T. Higuchi, M. Haruyama, K. Yamada, S. Onoda, T. Ohshima, W. Kada, O. Hanaizumi, A. Stacey, T. Teraji, S. Kono, J. Isoya and H. Kawarada, J. Phys. Chem. C 123, 3594 (2019).
  * [11] M. Kaviani, P. Deák, B. Aradi, T. Frauenheim, J.-P. Chou, and A. Gali, Nano Lett. 14, 4772 (2014).
  * [12] J.-P. Chou, A. Retzker, and A. Gali, Nano Lett. 17, 4, 2294 (2017).
  * [13] S. Li, J.-P. Chou, J. Wei, M. Sun, A. Hu, A. Gali, Carbon 145, 273 (2019).
  * [14] A. Gali, Nanophotonics, 8(11), 1907, (2019).
  * [15] C. Müller, X. Kong, J.-M. Cai, K. Melentijevic, A. Stacey, M. Markham, D. Twitchen, J. Isoya, S. Pezzagna, J. Meijer, J.F. Du, M.B. Plenio, B. Naydenov, L.P. McGuinness and F. Jelezko, Nat. Commun. 5, 4703 (2014).
  * [16] B. K. Ofori-Okai, S. Pezzagna,K. Chang, M. Loretz, R. Schirhagl, Y. Tao, B. A. Moores, K. Groot-Berning, J. Meijer, C. L. Degen, Phys. Rev. B 86, 081406(R) (2012).
  * [17] W. Körner, D. F. Urban and C. Elsässer, Phys. Rev. B 103, 085305 (2021).
  * [18] H. Holloway, K. C. Hass, M. A. Tamor, T. R. Anthony, and W. F. Banholzer, Phys. Rev. B 44, 7123 (1991).
  * [19] L. S. Ramsdell, Journal of Earth and Planetary Mat., 32, 64 (1947).
  * [20] A. Stacey, K. M. O’Donnell, J.-P. Chou, A. Schenk, A. Tadich, N. Dontschuk , J. Cervenka, C. Pakes, A.Gali, A. Hoffman, and S. Prawer, Adv. Mater. Interfaces 2, 1500079 (2015).
  * [21] G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * [22] G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * [23] P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * [24] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * [25] J. B. Cui, J. Ristein, and L. Ley, Phys. Rev. Lett. 81, 429 (1998).
  * [26] S. J. Sque, R. Jones and P. R. Briddon, Phys. Rev. B 73, 085313 (2006).
  * [27] A. K. Tiwari, J. P. Goss, P. R. Briddon, N. G. Wright, A. B. Horsfall, R. Jones, H. Pinto and M. J. Rayson, Phys. Rev. B 84, 245305 (2011).
  * [28] J.-P. Chou, and A. Gali, MRS Commun. 7, 551 (2017).
  * [29] D. Zhu, J. A. Bandy, S. Li, R. J. Hamers, Surf. Sci. 295, 650 (2016).
  * [30] A. Gali, M. Fyta and E. Kaxiras, Phys. Rev. B 77, 155206 (2008).
  * [31] S. Felton, A. M. Edmonds, M. E. Newton, P. M. Martineau, D. Fisher, D. J. Twitchen and J. M. Baker, Phys. Rev. B 79, 075203 (2009).
  * [32] T. Biktagirov, W. G. Schmidt and U. Gerstmann, Phys. Rev. Res. 2, 022024(R) (2020).
  * [33] R. C. Bansal, F. J. Vastola, and P. L. Walker, Carbon 10, 443 (1972).
  * [34] V. Y. Osipov, A. I. Shames, and A. Y. Vul, Phys. B Condens. Matter 404, 4522 (2009).
  * [35] Y. Romach, C. Müller, T. Unden, L. J. Rogers, T. Isoda, K. M. Itoh, M. Markham, A. Stacey, J. Meijer, S. Pezzagna, B. Naydenov, L. P. McGuinness, N. Bar-Gill, and F. Jelezko, Phys. Rev. Lett. 114, 017601 (2015).
  * [36] J. H. N. Loubser and J. A. van Wyk, Rep. Prog. Phys. 41, 1201 (1978).
  * [37] V. Ivády, T. Simon, J.R. Maze, I. A. Abrikosov and A. Gali, Phys. Rev. B 90, 235205 (2014).
  * [38] M. A. Bernstein; K. F. King; X. J. Zhou (2004). Handbook of MRI Pulse Sequences. San Diego: Elsevier Academic Press. p. 960. ISBN 0-12-092861-2.
  * [39] S. Kawai,H. Yamano,T. Sonoda, K. Kato, J. J. Buendia, T. Kageura, R. Fukuda, T. Okada, T. Tanii, T. Higuchi, M. Haruyama, K. Yamada, S. Onoda, T. Ohshima, W. Kada, O. Hanaizumi, A. Stacey,T. Teraji, S. Kono, J. Isoya, H. Kawarada, J. Phy. Chem. C,123,6 (2019).
  * [40] A. Stacey, K. O’Donnell, J. Chou, A. Schenk, A. Tadich, N. Dontschuk, J. Cervenka, K. Pakes, A. Gali, A. Hoffman, S. Prawer, Advanced Materials Interfaces 2, 10 (2015).