# PySINDy: A Python package for the Sparse Identification of Nonlinear
Dynamics from Data

Brian M. de Silva1∗, Kathleen Champion1∗, Markus Quade2,
Jean-Christophe Loiseau3, J. Nathan Kutz1, Steven L. Brunton4
1 Department of Applied Mathematics, University of Washington, Seattle, WA
98195, United States
2 Ambrosys GmbH, Potsdam, Germany
3 ParisTech, Paris, France
4 Department of Mechanical Engineering, University of Washington, Seattle, WA
98195, United States

###### Abstract

PySINDy is a Python package for the discovery of governing dynamical systems
models from data. In particular, PySINDy provides tools for applying the
sparse identification of nonlinear dynamics (SINDy) [1] approach to model
discovery. In this work we provide a brief description of the mathematical
underpinnings of SINDy, an overview and demonstration of the features
implemented in PySINDy (with code examples), practical advice for users, and a
list of potential extensions to PySINDy. Software is available at
https://github.com/dynamicslab/pysindy.

_Keywords–_ system identification, dynamical systems, symbolic regression,
open source, python

††∗ Corresponding authors<EMAIL_ADDRESS>and kpchamp@uw.edu).

##  1 Introduction

Scientists have long quantified empirical observations by developing
mathematical models that characterize the observations, have some measure of
interpretability, and are capable of making predictions. Dynamical systems
models, in particular, have been widely used to study, explain, and predict
behavior in a diversity of application areas, with examples ranging from
Newton’s laws of classical mechanics to the Michaelis-Menten kinetics for
modeling enzyme kinetics. While governing laws and equations have
traditionally been derived from first principles and expert knowledge, the
current growth of available measurement data and the resulting emphasis on
data-driven modeling motivates algorithmic and reproducible approaches for
automated model discovery.

A number of such approaches have been developed in recent years [2], including
linear methods [3, 4], dynamic mode decomposition (DMD) [5, 6] and Koopman
theory more generally [7, 8, 9, 10, 11, 12], nonlinear autoregressive
algorithms [13, 14], neural networks [15, 16, 17, 18, 19, 20, 21, 22, 23, 24],
Gaussian process regression [25, 26], operator inference and reduced-order
modeling [27, 28, 29], nonlinear Laplacian spectral analysis [30], diffusion
maps [31], genetic programming [32, 33, 34], and sparse regression [1, 35].
Maximizing the impact of these model discovery methods requires tools to make
them widely accessible to scientists across domains and at various levels of
mathematical expertise.

PySINDy is a Python package for the discovery of governing dynamical systems
models from data. In particular, PySINDy provides tools for applying the
sparse identification of nonlinear dynamics (SINDy) approach to model
discovery [1]. SINDy poses model discovery as a sparse regression problem,
where relevant terms in the dynamics are selected from a library of candidate
functions, many of them motivated by our deep historical knowledge of diverse
physics models. This approach results in interpretable models, and it has been
widely applied [36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50]
and extended [51, 37, 52, 53, 35, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63]
using different sparse optimization algorithms and library functions.

The PySINDy package is aimed at researchers and practitioners alike, enabling
anyone with access to measurement data to engage in scientific model
discovery. The package is designed to be accessible to inexperienced users,
adhering to scikit-learn standards, while also including customizable options
for more advanced users. A number of popular SINDy variants are implemented,
but PySINDy is also designed to enable further extensions for research and
experimentation.

##  2 Background

PySINDy provides an implementation of the SINDy method to discover governing
dynamical systems models of the form

$\frac{d}{dt}\mathbf{x}(t)=\mathbf{f}(\mathbf{x}(t)).$ (1)

Given data in the form of state measurements $\mathbf{x}(t)\in\mathbb{R}^{n}$,
SINDy identifies a model for the dynamics, given by the function $\mathbf{f}$,
which describes how the state of the system evolves in time. In particular,
SINDy sparsely approximates the dynamics in a library of candidate basis
functions
$\boldsymbol{\theta}(\mathbf{x})=[\theta_{1}(\mathbf{x}),\theta_{2}(\mathbf{x}),\dots,\theta_{\ell}(\mathbf{x})]$,
so that

$\mathbf{f}(\mathbf{x})\approx\sum_{k=1}^{\ell}\theta_{k}(\mathbf{x})\xi_{k}.$
(2)

The majority of coefficients $\xi_{k}$ are zero, and nonzero entries identify
active terms in the dynamics.

To pose SINDy as a regression problem, time-series measurements of
$\mathbf{x}$ and their time derivatives $\dot{\mathbf{x}}$ are arranged into
matrices

$\mathbf{X}=\left(\begin{array}[]{cccc}x_{1}(t_{1})&x_{2}(t_{1})&\cdots&x_{n}(t_{1})\\\
x_{1}(t_{2})&x_{2}(t_{2})&\cdots&x_{n}(t_{2})\\\
\vdots&\vdots&\ddots&\vdots\\\
x_{1}(t_{m})&x_{2}(t_{m})&\cdots&x_{n}(t_{m})\\\
\end{array}\right),\quad\dot{\mathbf{X}}=\left(\begin{array}[]{cccc}\dot{x}_{1}(t_{1})&\dot{x}_{2}(t_{1})&\cdots&\dot{x}_{n}(t_{1})\\\
\dot{x}_{1}(t_{2})&\dot{x}_{2}(t_{2})&\cdots&\dot{x}_{n}(t_{2})\\\
\vdots&\vdots&\ddots&\vdots\\\
\dot{x}_{1}(t_{m})&\dot{x}_{2}(t_{m})&\cdots&\dot{x}_{n}(t_{m})\\\
\end{array}\right).$

The derivatives can be approximated numerically or measured directly. The
library functions are evaluated on the data, resulting in
${\boldsymbol{\Theta}(\mathbf{X})=[\theta_{1}(\mathbf{X}),\theta_{2}(\mathbf{X}),\dots,\theta_{\ell}(\mathbf{X})]}$.
Sparse regression is then performed to approximately solve

$\dot{\mathbf{X}}\approx\boldsymbol{\Theta}(\mathbf{X})\boldsymbol{\Xi},$ (3)

where $\boldsymbol{\Xi}$ is a set of coefficients that determines the active
terms in $\mathbf{f}$. While the original SINDy formulation solves (3) using a
sequentially thresholded least squares algorithm [1, 64], this problem can be
solved using any sparse regression algorithm, such as lasso [65], sparse
relaxed regularized regression (SR3) [66, 67], stepwise sparse regression
(SSR) [58], or Bayesian methods [68, 69, 70].

SINDy has been widely applied for model identification in applications such as
chemical reaction dynamics [40], nonlinear optics [36], fluid dynamics [37,
39, 41, 42, 47] and turbulence modeling [48, 50], plasma convection [38],
numerical algorithms [45], and structural modeling [46], among others [43, 44,
49]. It has also been extended to handle more complex modeling scenarios such
as partial differential equations [52, 35], systems with inputs or control
[59], systems with implicit dynamics [51], hybrid systems [60], to enforce
physical constraints [37], to incorporate information theory [53], to identify
models from corrupt or limited data [54, 56] and ensembles of initial
conditions [57], and extending the formulation to include integral terms [55,
63], tensor representations [61, 62], and stochastic forcing [58]. However,
there is not a definitive standard implementation or package for SINDy.
Versions of SINDy have been implemented within larger projects such as
sparsereg [71], but no specific implementation has emerged as the most widely
adopted and most versions implement a limited set of features. Researchers
have thus typically written their own implementations, resulting in duplicated
effort and a lack of standardization. This not only makes it more difficult to
apply SINDy to scientific data sets, but also makes it more challenging to
benchmark extensions of the method and makes such extensions less accessible
to end users. This motivates the creation of a dedicated package for SINDy.
The PySINDy package provides a central codebase where many of the basic SINDy
features are implemented, allowing for easy use and standardization. In
addition, it is straightforward for users to extend PySINDy so that new
developments are available to the wider community.

##  3 Features

The core object in the PySINDy package is the SINDy model class, which is
implemented as a scikit-learn estimator. This design choice was made to ensure
that the package is simple to use for a wide user base, as many potential
users will be familiar with scikit-learn. It also expresses the SINDy model
object at the appropriate level of abstraction so that users can embed it into
more sophisticated pipelines in scikit-learn, such as for parameter tuning and
model selection.

Our PySINDy implementation involves three major steps, resulting in three
modeling decisions:

  1. 1.

The numerical differentiation scheme used to approximate $\dot{\mathbf{X}}$
from $\mathbf{X}$;

  2. 2.

The candidate functions constituting the feature library
$\boldsymbol{\Theta}$;

  3. 3.

The sparse regression algorithm that is applied to solve (3) to find
$\boldsymbol{\Xi}$.

$\dot{\mathbf{X}}$ $\approx$ $\boldsymbol{\Theta}(\mathbf{X})$
$\boldsymbol{\Xi}$ pysindy.feature_library pysindy.differentiation
pysindy.optimizers Numerical differentiationForm feature librarySparse
regression Figure 1: Correspondence between the sparse regression problem
solved by SINDy and the submodules of PySINDy.

The core SINDy object was designed to incorporate these three components in as
modular a manner as possible, having one attribute corresponding to each:
SINDy.differentiation_method for numerical differentiation,
SINDy.feature_library for the formation of the candidate function library, and
SINDy.optimizer for the sparse regressor. PySINDy provides standard options
for each step, while making it easy to replace any of these steps with more
sophisticated “third-party” algorithms. In particular, at the time of writing,
we have implemented the following methods:

  * •

Numerical differentiation (for computing $\dot{\mathbf{X}}$ from $\mathbf{X}$)

    * –

Finite difference: FiniteDifference

    * –

Smoothed finite difference: SmoothedFiniteDifference

  * •

Feature library (for constructing $\boldsymbol{\Theta}$)

    * –

Multivariate polynomials: PolynomialLibrary

    * –

Fourier modes (i.e. trigonometric functions): FourierLibrary

    * –

Custom library (defined by user-supplied functions): CustomLibrary

    * –

Identity library (in case users want to compute $\boldsymbol{\Theta}$
themselves): IdentityLibrary

  * •

Optimizer (for performing sparse regression)

    * –

Sequentially thresholded least-squares [1, 64]: STLSQ

    * –

Sparse relaxed regularized regression (SR3) [66]: SR3

##  4 Examples

The PySINDy GitHub page111https://github.com/dynamicslab/pysindy includes
tutorials in the form of Jupyter notebooks. These tutorials demonstrate the
usage of various features of the package and reproduce the examples from the
original SINDy paper [1]. Throughout this section we will use the Lorenz
equations (4) as the dynamical system to illustrate the PySINDy package:

$\left\\{\begin{aligned} \dot{x}&=-10x-10y\\\ \dot{y}&=x(28-z)-y\\\
\dot{z}&=xy-\frac{8}{3}z\end{aligned}\right.$ (4)

In Python, the right-hand side of (4) can be expressed as follows:

⬇

def lorenz(x, t):

return [

10 ∗ (x[1] - x[0]),

x[0] ∗ (28 - x[2]) - x[1],

x[0] ∗ x[1] - (8 / 3) ∗ x[2]

]

To construct training data to feed into a SINDy model, we integrate (4) with:

⬇

import numpy as np

from scipy.integrate import odeint

dt = 0.002

t = np.arange(0, 10, dt)

x0 = [-8, 8, 27]

X = odeint(lorenz, x0, t)

We plot X in Figure 2. It is important to note that each column of X
corresponds to a variable and each row to a point in time. All PySINDy objects
that handle data assume the data is structured this way.

Figure 2: Measurement data simulated using the Lorenz equations (4).

###  4.1 Basic usage

The pysindy package is built around the SINDy class, which encapsulates all
the steps necessary to learn a dynamical system with SINDy. To create a SINDy
object, fit it to the data (i.e. to infer a dynamical system from the data),
and print the resulting model, we invoke the SINDy constructor, the fit
method, and custom print functions

⬇

model = ps.SINDy()

model.fit(X, t=dt)

model.print()

which generates the following output

      x0’ = -9.999 x0 + 9.999 x1
      x1’ = 27.992 x0 + -0.999 x1 + -1.000 x0 x2
      x2’ = -2.666 x2 + 1.000 x0 x1

Once the SINDy object has been fit we can feed in new data and use the learned
model to predict the derivatives for each measurement (recall that
measurements correspond to rows).

⬇

t_test = np.arange(0, 15, dt)

x0_test = np.array([8, 7, 15])

X_test = odeint(lorenz, x0_test, t_test)

X_dot_test_computed = model.differentiate(X_test, t=dt)

X_dot_test_predicted = model.predict(X_test)

The call model.differentiate(X_test, t=dt) applies the numerical
differentiation method in the SINDy model to X_test with time steps of length
dt. In Figure 3 we plot each dimension of X_dot_test_computed and
X_dot_test_predicted.

Figure 3: Derivatives of variables from the Lorenz equation via numerical
differentiation and using a learned SINDy model.

Rather than predicting derivatives, we will often be interested in using our
model to evolve initial conditions forward in time using the learned model.
The simulate function does just that.

⬇

X_test_sim = model.simulate(x0_test, t_test)

Figure 4 shows the simulated trajectory plotted against the true trajectory
X_test. The trajectories agree initially, but they eventually diverge due to
the chaotic nature of the Lorenz equations.

Figure 4: Two trajectories starting at the same position evolved forward in
time with the exact Lorenz equations (black, solid) and the learned SINDy
model (red, dashed).

###  4.2 Custom features

Thus far we have relied on the default options of the SINDy object, but
PySINDy comes equipped with multiple alternative built-in methods for
differentiation, library building, and optimization. These options are
selected by passing corresponding PySINDy objects to the SINDy constructor via
the differentiation_method, feature_library, and optimizer arguments,
respectively. Parameters for the differentiation, library, and optimization
algorithms are supplied to the corresponding objects’ constructors rather than
directly to the SINDy object. We demonstrate the syntax with the following
example.

⬇

differentiation_method = ps.FiniteDifference(order=1)

feature_library = ps.PolynomialLibrary(degree=3, include_bias=False)

optimizer = ps.SR3(threshold=0.1, nu=1, tol=1e-6)

model = ps.SINDy(

differentiation_method=differentiation_method,

feature_library=feature_library,

optimizer=optimizer,

feature_names=["x", "y", "z"]

)

model.fit(X, t=dt)

model.print()

which prints

      x’ = -10.021 x + 9.993 y
      y’ = 28.431 x + -1.212 y + -1.008 x z
      z’ = -2.675 z + 1.000 x y.

A number of other built-in options are available. The official
documentation222https://pysindy.readthedocs.io/en/latest/index.html and
examples333https://github.com/dynamicslab/pysindy/tree/master/example provide
an exhaustive list.

##  5 Practical tips

In this section we provide pragmatic advice for using PySINDy effectively. We
discuss potential pitfalls and strategies for overcoming them. We also specify
how to incorporate custom methods not implemented natively in PySINDy, where
applicable. The information presented here is derived from a combination of
experience and theoretical considerations.

###  5.1 Numerical differentiation

Numerical differentiation is one of the core components of the SINDy method.
Derivatives of measurement variables provide the targets for the sparse
regression problem (3). If care is not taken in computing these derivatives,
the quality of the learned model is likely to suffer.

By default, a second order finite difference method is used to differentiate
input data. Finite difference methods tend to amplify noise in data. If the
data are smooth (at least twice differentiable), then finite difference
methods give accurate derivative approximations. When the data are noisy, they
give derivative estimates with more noise than the original data. Figure 5
visualizes the impact of noise on numerical derivatives. Note that even a
small amount of noise in the data can produce noticeable degradation in the
quality of the numerical derivative.

Figure 5: A toy example illustrating the effect of noise on derivatives
computed with a second order finite difference method. Left: The data to be
differentiated; $y=\sin(x)$ with and without a small amount of additive noise
(normally distributed with mean 0 and standard deviation $0.01$). Right:
Derivatives of the data; the exact derivative $\cos(x)$ (blue), the finite
difference derivative of the exact data (black, dashed), and the finite
difference derivative of the noisy data.

One way to mitigate the effects of noise is to smooth the measurements before
computing derivatives. The SmoothedFiniteDifference method can be used for
this purpose. A numerical differentiation scheme with total variation
regularization has also been proposed [72] and recommended for use in SINDy
[1].

Users wishing to employ their own numerical differentiation schemes have two
ways of doing so. Derivatives of input measurements can be computed externally
with the method of choice and then passed directly into the SINDy.fit method
via the x_dot keyword argument. Alternatively, users can implement their own
differentiation methods and pass them into the SINDy constructor using the
differentiation_method argument. In this case, the supplied class need only
have implemented a __call__ method taking two arguments, x and t.

###  5.2 Library selection

The SINDy method assumes dynamics can be represented as a sparse linear
combination of library functions. If this assumption is violated, the method
is likely to exhibit poor performance. This issue tends to manifest itself as
numerous library terms being active, often with weights of vastly different
magnitudes, still resulting in poor model error.

Typically, prior knowledge of the system of interest and its dynamics should
be used to make a judicious choice of basis functions. When such information
is unavailable, the default class of library functions, polynomials, are a
good place to start, as smooth functions have rapidly converging Taylor
series. Brunton et al. [1] showed that, equipped with a polynomial library,
SINDy can recover the first few terms of the (zero-centered) Taylor series of
the true right-hand side function $\mathbf{f}(x)$. If one has reason to
believe the dynamics can be sparsely represented in terms of Chebyshev
polynomials rather than monomials, then the library should include Chebyshev
polynomials.

PySINDy includes the CustomLibrary and IdentityLibrary objects to allow for
flexibility in the library functions. When the desired library consists of a
set of functions that should be applied to each measurement variable (or pair,
triplet, etc. of measurement variables) in turn, the CustomLibrary class
should be used. The IdentityLibrary class is the most customizable, but
transfers the work of computing library functions over to the user. It expects
that all the features one wishes to include in the library have already been
computed and are present in X before SINDy.fit is called, as it simply applies
the identity map to each variable that is passed to it. It is best suited for
situations in which one has very specific instructions for how to apply
library functions (e.g. if some of the functions should be applied to only
some of the input variables).

As terms are added to the library, the underlying sparse regression problem
becomes increasingly ill-conditioned. Therefore it is recommended to start
with a small library whose size is gradually expanded until the desired level
of performance is achieved. For example, a user may wish to start with a
library of linear terms and then add quadratic and cubic terms as necessary to
improve model performance. For the best results, the strength of
regularization applied should be increased in proportion to the size of the
library to account for the worsening condition number of the resulting linear
system.

Users may also choose to implement library classes tailored to their
applications. To do so one should have the new class inherit from our
BaseFeatureLibrary class. See the documentation for guidance on which
functions the new class is expected to implement.

###  5.3 Optimization

PySINDy uses various optimizers to solve the sparse regression problem. For a
fixed differentiation method, set of inputs, and candidate library, there is
still some variance in the dynamical system identified by SINDY, depending on
which optimizer is employed.

The default optimizer in PySINDy is the sequentially-thresholded least-squares
algorithm (STLSQ). In addition to being the method originally proposed for use
with SINDy, it involves a single, easily interpretable hyperparameter, and it
exhibits good performance across a variety of problems.

The sparse relaxed regularized regression (SR3) [66, 67] algorithm can be used
when the results of STLSQ are unsatisfactory. It involves a few more
hyperparameters that can be tuned for improved accuracy. In particular, the
thresholder parameter controls the type of regularization that is applied. For
optimal results, one may find it useful to experiment with $L^{0}$, $L^{1}$,
and clipped absolute deviation (CAD) regularization. The other hyperparameters
can be tuned with cross-validation.

Custom or third party sparse regression methods are also supported. Simply
instantiate an instance of the custom object and pass it to the SINDy
constructor using the optimizer keyword. Our implementation is compatible with
any of the linear models from Scikit-learn (e.g. RidgeRegression, Lasso, and
ElasticNet). See the documentation for a list of methods and attributes a
custom optimizer is expected to implement. There you will also find an example
where the Scikit-learn Lasso object is used to perform sparse regression.

###  5.4 Regularization

Regularization, in this context, is a technique for improving the conditioning
of ill-posed problems. Without regularization, one often obtains highly
unstable results, with learned parameter values differing substantially for
slightly different inputs. SINDy seeks weights that express dynamics as a
sparse linear combination of library functions. When the columns of the
measurement data or the library are statistically correlated, which is likely
for large libraries, the SINDy inverse problem can quickly become ill-posed.
Though the sparsity constraint is a type of regularization itself, for many
problems another form of regularization is needed for SINDy to learn a robust
dynamical model.

In some cases regularization can be interpreted as enforcing a prior
distribution on the model parameters [73]. Applying strong regularization
biases the learned weights away from the values that would allow them to best
fit the data and toward the values preferred by the prior distribution (e.g.
$L^{2}$ regularization corresponds to a Gaussian prior). Therefore once a
sparse set of nonzero coefficients is discovered, our methods apply an extra
“unbiasing” step where unregularized least-squares is used to find the values
of the identified nonzero coefficients. All of our built-in methods use
regularization by default.

Some general best practices regarding regularization follow. Most problems
will benefit from some amount of regularization. Regularization strength
should be increased as the size of the candidate right-hand side library
grows. If warnings about ill-conditioned matrices are generated when SINDy.fit
is called, more regularization may help. We also recommend setting unbias to
True when invoking the SINDy.fit method, especially when large amounts of
regularization are being applied. Cross-validation can be used to select
appropriate regularization parameters for a given problem.

##  6 Extensions

In this section we list potential extensions and enhancements to our SINDy
implementation. We provide references for the improvements that are inspired
by previously conducted research and the rationale behind the other potential
changes.

  * •

Partial differential equations (PDEs): While dynamical systems given by
ordinary differential equations (ODEs) provide a flexible approach to modeling
physical systems, many systems are inherently described by partial
differential equations (PDEs), which are not immediately discoverable using
PySINDy. Multiple approaches for the data-driven discovery of PDEs have been
proposed [35, 52] as extensions to the SINDy method, and these may be readily
included within the PySINDy framework.

  * •

Identifying coordinates and latent variables: Many complex systems, such as
fluid flows, are high-dimensional, yet exhibit low-dimensional patterns that
may be exploited for modeling [74, 75, 76]. Identifying effective coordinate
systems on which to build models is an important aspect of data-driven
discovery. Recently, SINDy has been embedded into an autoencoder framework
[23] to simultaneously identify effective coordinates and sparse dynamics.
Similarly, for many systems, it is impossible to measure the full state of the
system, so that there are latent variables. Time-delay coordinates have been
useful for identifying sparse models from limited measurements [12]. Both of
these are candidates for future extensions.

  * •

Constraints: SINDy has been extended to enforce physical constraints during
the sparse regression step [37]. When working with physical systems with known
conserved quantities, such as the conservation of energy in incompressible
fluids [37], such a method allows one to automatically incorporate this prior
information into the model discovery process.

  * •

Integral formulation: We previously discussed how measurement data with too
much noise can disrupt the model discovery process, and we offered smoothing
as one possible solution. Another is to work with an integral version of (1),
as proposed by Schaeffer and McCalla [55] and extended to PDEs by Reinbold,
Gurevich, and Grigoriev [63]. Where numerical differentiation tends to amplify
noise, numerical integration tends to smooth it out. This formulation has been
shown to improve the robustness of SINDy to noise.

  * •

Ensembles: Ensembles are a proven class of methods in machine learning for
variance reduction (improved model generalizability) at the cost of extra
computation. Rather than training a single model, multiple high-variance
models are trained and their predictions are averaged together. We think that
ideas from ensemble learning could be adapted to improve the performance of
SINDy models. Recent work has hinted at possible approaches [77].

  * •

Extended libraries: Choosing the appropriate basis in which to represent
dynamics is of critical importance for the successful application of SINDy.
Although we currently provide methods allowing users the flexibility to input
their own library functions, we aim to make the library construction process
even easier by providing a common suite of tools for the creation and
combination of sets of candidate functions. More basis functions could be
supported natively, such as nonautonomous terms (those depending explicitly on
the dependent variable, time). Taking this idea a step further, specific
variables (columns of $\mathbf{X}$) which should not appear on the left-hand
side of (3) could be identified by the user. This would enable SINDy to
include inputs and control variables [59]. Operations acting on one or more
libraries could also be implemented. For example, combining libraries via
union, intersection, composition, or tensor product could enable the
expression of complicated nonlinear dynamics. The ability to apply libraries
to only subsets of state variables could help cut down on the computational
cost and improve the conditioning of the sparse regression problem solved
within SINDy.

##  7 Acknowledgments

PySINDy is a fork of sparsereg [71]. SLB acknowledges funding support from the
Army Research Office (ARO W911NF-19-1-0045) and the Air Force Office of
Scientific Research (AFOSR FA9550- 18-1-0200). JNK acknowledges support from
the Air Force Office of Scientific Research (AFOSR FA9550-17-1-0329). This
material is based upon work supported by the National Science Foundation
Graduate Research Fellowship under Grant Number DGE- 1256082.

## References

  * [1] S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Discovering governing equations from data by sparse identification of nonlinear dynamical systems,” _Proceedings of the National Academy of Sciences_ , vol. 113, no. 15, pp. 3932–3937, 2016.
  * [2] S. L. Brunton and J. N. Kutz, _Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control_. Cambridge University Press, 2019.
  * [3] O. Nelles, _Nonlinear system identification: from classical approaches to neural networks and fuzzy models_. Springer, 2013.
  * [4] L. Ljung, “Perspectives on system identification,” _Annual Reviews in Control_ , vol. 34, no. 1, pp. 1–12, 2010.
  * [5] P. J. Schmid, “Dynamic mode decomposition of numerical and experimental data,” _Journal of fluid mechanics_ , vol. 656, pp. 5–28, 2010.
  * [6] J. N. Kutz, S. L. Brunton, B. W. Brunton, and J. L. Proctor, _Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems_. SIAM, 2016.
  * [7] M. Budišić, R. Mohr, and I. Mezić, “Applied Koopmanism a),” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 22, no. 4, p. 047510, 2012.
  * [8] I. Mezic, “Analysis of fluid flows via spectral properties of the Koopman operator,” _Annual Review of Fluid Mechanics_ , vol. 45, pp. 357–378, 2013\.
  * [9] M. O. Williams, I. G. Kevrekidis, and C. W. Rowley, “A data-driven approximation of the Koopman operator: extending dynamic mode decomposition,” _Journal of Nonlinear Science_ , vol. 6, pp. 1307–1346, 2015\.
  * [10] S. Klus, F. Nüske, P. Koltai, H. Wu, I. Kevrekidis, C. Schütte, and F. Noé, “Data-driven model reduction and transfer operator approximation,” _Journal of Nonlinear Science_ , 2018.
  * [11] Q. Li, F. Dietrich, E. M. Bollt, and I. G. Kevrekidis, “Extended dynamic mode decomposition with dictionary learning: A data-driven adaptive spectral decomposition of the koopman operator,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 27, no. 10, p. 103111, 2017.
  * [12] S. L. Brunton, B. W. Brunton, J. L. Proctor, E. Kaiser, and J. N. Kutz, “Chaos as an intermittently forced linear system,” _Nature Communications_ , vol. 8, no. 19, pp. 1–9, 2017.
  * [13] H. Akaike, “Fitting autoregressive models for prediction,” _Ann Inst Stat Math_ , vol. 21, no. 1, pp. 243–247, 1969.
  * [14] S. A. Billings, _Nonlinear system identification: NARMAX methods in the time, frequency, and spatio-temporal domains_. John Wiley & Sons, 2013.
  * [15] Z. Long, Y. Lu, X. Ma, and B. Dong, “Pde-net: Learning pdes from data,” _arXiv preprint arXiv:1710.09668_ , 2017.
  * [16] L. Yang, D. Zhang, and G. E. Karniadakis, “Physics-informed generative adversarial networks for stochastic differential equations,” _arXiv preprint arXiv:1811.02033_ , 2018.
  * [17] C. Wehmeyer and F. Noé, “Time-lagged autoencoders: Deep learning of slow collective variables for molecular kinetics,” _The Journal of Chemical Physics_ , vol. 148, no. 241703, pp. 1–9, 2018.
  * [18] A. Mardt, L. Pasquali, H. Wu, and F. Noé, “VAMPnets: Deep learning of molecular kinetics,” _Nature Communications_ , vol. 9, no. 5, 2018.
  * [19] P. R. Vlachas, W. Byeon, Z. Y. Wan, T. P. Sapsis, and P. Koumoutsakos, “Data-driven forecasting of high-dimensional chaotic systems with long short-term memory networks,” _Proc. R. Soc. A_ , vol. 474, no. 2213, p. 20170844, 2018.
  * [20] J. Pathak, B. Hunt, M. Girvan, Z. Lu, and E. Ott, “Model-free prediction of large spatiotemporally chaotic systems from data: a reservoir computing approach,” _Physical review letters_ , vol. 120, no. 2, p. 024102, 2018.
  * [21] L. Lu, X. Meng, Z. Mao, and G. E. Karniadakis, “Deepxde: A deep learning library for solving differential equations,” _arXiv preprint arXiv:1907.04502_ , 2019.
  * [22] M. Raissi, P. Perdikaris, and G. Karniadakis, “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,” _Journal of Computational Physics_ , vol. 378, pp. 686–707, 2019.
  * [23] K. Champion, B. Lusch, J. N. Kutz, and S. L. Brunton, “Data-driven discovery of coordinates and governing equations,” _Proceedings of the National Academy of Sciences_ , vol. 116, no. 45, pp. 22 445–22 451, 2019.
  * [24] M. Raissi, A. Yazdani, and G. E. Karniadakis, “Hidden fluid mechanics: Learning velocity and pressure fields from flow visualizations,” _Science_ , vol. 367, no. 6481, pp. 1026–1030, 2020.
  * [25] M. Raissi and G. E. Karniadakis, “Machine learning of linear differential equations using gaussian processes,” _arXiv preprint arXiv:1701.02440_ , 2017\.
  * [26] ——, “Hidden physics models: Machine learning of nonlinear partial differential equations,” _Journal of Computational Physics_ , vol. 357, pp. 125–141, 2018.
  * [27] P. Benner, S. Gugercin, and K. Willcox, “A survey of projection-based model reduction methods for parametric dynamical systems,” _SIAM review_ , vol. 57, no. 4, pp. 483–531, 2015.
  * [28] B. Peherstorfer and K. Willcox, “Data-driven operator inference for nonintrusive projection-based model reduction,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 306, pp. 196–215, 2016.
  * [29] E. Qian, B. Kramer, B. Peherstorfer, and K. Willcox, “Lift & learn: Physics-informed machine learning for large-scale nonlinear dynamical systems,” _Physica D: Nonlinear Phenomena_ , vol. 406, p. 132401, 2020.
  * [30] D. Giannakis and A. J. Majda, “Nonlinear laplacian spectral analysis for time series with intermittency and low-frequency variability,” _Proceedings of the National Academy of Sciences_ , vol. 109, no. 7, pp. 2222–2227, 2012.
  * [31] O. Yair, R. Talmon, R. R. Coifman, and I. G. Kevrekidis, “Reconstruction of normal forms by learning informed observation geometries from data,” _Proceedings of the National Academy of Sciences_ , p. 201620045, 2017.
  * [32] J. Bongard and H. Lipson, “Automated reverse engineering of nonlinear dynamical systems,” _Proceedings of the National Academy of Sciences_ , vol. 104, no. 24, pp. 9943–9948, 2007.
  * [33] M. Schmidt and H. Lipson, “Distilling free-form natural laws from experimental data,” _Science_ , vol. 324, no. 5923, pp. 81–85, 2009.
  * [34] B. C. Daniels and I. Nemenman, “Automated adaptive inference of phenomenological dynamical models,” _Nature communications_ , vol. 6, 2015\.
  * [35] S. H. Rudy, S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Data-driven discovery of partial differential equations,” _Science Advances_ , vol. 3, no. e1602614, 2017.
  * [36] M. Sorokina, S. Sygletos, and S. Turitsyn, “Sparse identification for nonlinear optical communication systems: SINO method,” _Optics express_ , vol. 24, no. 26, pp. 30 433–30 443, 2016.
  * [37] J.-C. Loiseau and S. L. Brunton, “Constrained sparse Galerkin regression,” _Journal of Fluid Mechanics_ , vol. 838, pp. 42–67, 2018.
  * [38] M. Dam, M. Brøns, J. Juul Rasmussen, V. Naulin, and J. S. Hesthaven, “Sparse identification of a predator-prey system from simulation data of a convection model,” _Physics of Plasmas_ , vol. 24, no. 2, p. 022310, 2017\.
  * [39] J.-C. Loiseau, B. R. Noack, and S. L. Brunton, “Sparse reduced-order modeling: sensor-based dynamics to full-state estimation,” _Journal of Fluid Mechanics_ , vol. 844, pp. 459–490, 2018.
  * [40] M. Hoffmann, C. Fröhner, and F. Noé, “Reactive SINDy: Discovering governing reactions from concentration data,” _Journal of Chemical Physics_ , vol. 150, no. 025101, 2019.
  * [41] J.-C. Loiseau, “Data-driven modeling of the chaotic thermal convection in an annular thermosyphon,” _arXiv preprint arXiv:1911.07920_ , 2019.
  * [42] Y. El Sayed M, R. Semaan, and R. Radespiel, “Sparse modeling of the lift gains of a high-lift configuration with periodic coanda blowing,” in _2018 AIAA Aerospace Sciences Meeting_ , 2018, p. 1054.
  * [43] A. Narasingam and J. S.-I. Kwon, “Data-driven identification of interpretable reduced-order models using sparse regression,” _Computers & Chemical Engineering_, vol. 119, pp. 101–111, 2018.
  * [44] B. de Silva, D. M. Higdon, S. L. Brunton, and J. N. Kutz, “Discovery of physics from data: Universal laws and discrepancies,” _arXiv preprint arXiv:1906.07906_ , 2019.
  * [45] S. Thaler, L. Paehler, and N. A. Adams, “Sparse identification of truncation errors,” _Journal of Computational Physics_ , vol. 397, p. 108851, 2019.
  * [46] Z. Lai and S. Nagarajaiah, “Sparse structural system identification method for nonlinear dynamic systems with hysteresis/inelastic behavior,” _Mechanical Systems and Signal Processing_ , vol. 117, pp. 813–842, 2019\.
  * [47] N. Deng, B. R. Noack, M. Morzynski, and L. R. Pastur, “Low-order model for successive bifurcations of the fluidic pinball,” _Journal of fluid mechanics_ , vol. 884, no. A37, 2020.
  * [48] M. Schmelzer, R. P. Dwight, and P. Cinnella, “Discovery of algebraic reynolds-stress models using sparse symbolic regression,” _Flow, Turbulence and Combustion_ , vol. 104, no. 2, pp. 579–603, 2020.
  * [49] S. Pan, N. Arnold-Medabalimi, and K. Duraisamy, “Sparsity-promoting algorithms for the discovery of informative koopman invariant subspaces,” _arXiv preprint arXiv:2002.10637_ , 2020.
  * [50] S. Beetham and J. Capecelatro, “Formulating turbulence closures using sparse regression with embedded form invariance,” _arXiv preprint arXiv:2003.12884_ , 2020.
  * [51] N. M. Mangan, S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Inferring biological networks by sparse identification of nonlinear dynamics,” _IEEE Transactions on Molecular, Biological, and Multi-Scale Communications_ , vol. 2, no. 1, pp. 52–63, 2016.
  * [52] H. Schaeffer, “Learning partial differential equations via data discovery and sparse optimization,” in _Proc. R. Soc. A_ , vol. 473, no. 2197. The Royal Society, 2017, p. 20160446.
  * [53] N. M. Mangan, J. N. Kutz, S. L. Brunton, and J. L. Proctor, “Model selection for dynamical systems via sparse regression and information criteria,” _Proceedings of the Royal Society A_ , vol. 473, no. 2204, pp. 1–16, 2017\.
  * [54] G. Tran and R. Ward, “Exact recovery of chaotic systems from highly corrupted data,” _Multiscale Modeling & Simulation_, vol. 15, no. 3, pp. 1108–1129, 2017.
  * [55] H. Schaeffer and S. G. McCalla, “Sparse model selection via integral terms,” _Physical Review E_ , vol. 96, no. 2, p. 023302, 2017.
  * [56] H. Schaeffer, G. Tran, and R. Ward, “Extracting sparse high-dimensional dynamics from limited data,” _SIAM Journal on Applied Mathematics_ , vol. 78, no. 6, pp. 3279–3295, 2018.
  * [57] K. Wu and D. Xiu, “Numerical aspects for approximating governing equations using data,” _arXiv preprint arXiv:1809.09170_ , 2018.
  * [58] L. Boninsegna, F. Nüske, and C. Clementi, “Sparse learning of stochastic dynamical equations,” _The Journal of chemical physics_ , vol. 148, no. 24, p. 241723, 2018.
  * [59] E. Kaiser, J. N. Kutz, and S. L. Brunton, “Sparse identification of nonlinear dynamics for model predictive control in the low-data limit,” _Proceedings of the Royal Society of London A_ , vol. 474, no. 2219, 2018\.
  * [60] N. M. Mangan, T. Askham, S. L. Brunton, J. N. Kutz, and J. L. Proctor, “Model selection for hybrid dynamical systems via sparse regression,” _Proceedings of the Royal Society A_ , vol. 475, no. 2223, p. 20180534, 2019\.
  * [61] P. Gelß, S. Klus, J. Eisert, and C. Schütte, “Multidimensional approximation of nonlinear dynamical systems,” _Journal of Computational and Nonlinear Dynamics_ , vol. 14, no. 6, 2019.
  * [62] A. Goeßmann, M. Götte, I. Roth, R. Sweke, G. Kutyniok, and J. Eisert, “Tensor network approaches for learning non-linear dynamical laws,” _arXiv preprint arXiv:2002.12388_ , 2020.
  * [63] P. A. Reinbold, D. R. Gurevich, and R. O. Grigoriev, “Using noisy or incomplete data to discover models of spatiotemporal dynamics,” _Physical Review E_ , vol. 101, no. 1, p. 010203, 2020.
  * [64] L. Zhang and H. Schaeffer, “On the convergence of the sindy algorithm,” _Multiscale Modeling & Simulation_, vol. 17, no. 3, pp. 948–972, 2019.
  * [65] R. Tibshirani, “Regression shrinkage and selection via the lasso,” _Journal of the Royal Statistical Society. Series B (Methodological)_ , vol. 58, no. 1, pp. 267–288, 1996. [Online]. Available: http://www.jstor.org/stable/2346178
  * [66] P. Zheng, T. Askham, S. L. Brunton, J. N. Kutz, and A. Y. Aravkin, “A unified framework for sparse relaxed regularized regression: Sr3,” _IEEE Access_ , vol. 7, pp. 1404–1423, 2018.
  * [67] K. Champion, P. Zheng, A. Y. Aravkin, S. L. Brunton, and J. N. Kutz, “A unified sparse optimization framework to learn parsimonious physics-informed models from data,” _arXiv preprint arXiv:1906.10612_ , 2019.
  * [68] S. Zhang and G. Lin, “Robust data-driven discovery of governing physical laws with error bars,” _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , vol. 474, no. 2217, p. 20180305, 2018.
  * [69] W. Pan, Y. Yuan, J. Gonçalves, and G. Stan, “A sparse Bayesian approach to the identification of nonlinear state-space systems,” _IEEE Transactions on Automatic Control_ , vol. 61, no. 1, pp. 182–187, January 2016\.
  * [70] R. K. Niven, A. Mohammad-Djafari, L. Cordier, M. Abel, and M. Quade, “Bayesian identification of dynamical systems,” _Multidisciplinary Digital Publishing Institute Proceedings_ , vol. 33, no. 1, p. 33, 2020.
  * [71] M. Quade, “sparsereg - collection of modern sparse regression algorithms,” Feb. 2018. [Online]. Available: https://github.com/ohjeah/sparsereg
  * [72] R. Chartrand, “Numerical differentiation of noisy, nonsmooth data,” _ISRN Applied Mathematics_ , vol. 2011, 2011.
  * [73] C. M. Bishop, _Pattern recognition and machine learning_. springer, 2006.
  * [74] K. Taira, S. L. Brunton, S. Dawson, C. W. Rowley, T. Colonius, B. J. McKeon, O. T. Schmidt, S. Gordeyev, V. Theofilis, and L. S. Ukeiley, “Modal analysis of fluid flows: An overview,” _AIAA Journal_ , vol. 55, no. 12, pp. 4013–4041, 2017.
  * [75] K. Taira, M. S. Hemati, S. L. Brunton, Y. Sun, K. Duraisamy, S. Bagheri, S. Dawson, and C.-A. Yeh, “Modal analysis of fluid flows: Applications and outlook,” _AIAA Journal_ , vol. 58, no. 3, pp. 998–1022, 2020.
  * [76] S. L. Brunton, B. R. Noack, and P. Koumoutsakos, “Machine learning for fluid mechanics,” _Annual Review of Fluid Mechanics_ , vol. 52, pp. 477–508, 2020\.
  * [77] P. Sachdeva, J. Livezey, A. Tritt, and K. Bouchard, “Pyuoi: The union of intersections framework in python,” _Journal of Open Source Software_ , vol. 4, no. 44, p. 1799, 2019.


††thanks: Electronic address<EMAIL_ADDRESS>Corresponding author

#  _Ab initio_ study of (100) diamond surface spins

Jyh-Pin Chou  Department of Physics, National Changhua University of
Education, No.1, Jin-De Road, Changhua City, Changhua, Taiwan  Nathalie P. de
Leon  Department of Electrical and Computer Engineering, Princeton University,
Engineering Quadrangle, Olden Street, Princeton, NJ 08544, United States  Adam
Gali  Institute for Solid State Physics and Optics, Wigner Research Centre for
Physics, Budapest, POB 49, H-1525, Hungary  Department of Atomic Physics,
Budapest University of Technology and Economics, Budafoki út 8, H-1111,
Budapest, Hungary

###### Abstract

Unpaired electronic spins at diamond surfaces are ubiquitous and can lead to
decoherence of shallow nitrogen vacancy centers. They have been observed in
several studies to date, but their exact chemical nature is still unknown. We
propose a simple and reliable model to explain the existence and chemical
stability of surface spins on the (100) diamond surface by means of density
functional theory simulation. We find that the (111) facet, which is naturally
generated at a step edge of (100) crystalline diamond surface, can
topographically protect a spinful defect. Our study reveals a mechanism for
annihilation of these surface spins upon annealing. We also demonstrate that
the Fermi-contact term in the hyperfine coupling is not negligible between the
surface spins and the surrounding nuclear spins, and thus ab initio simulation
can be used to devise a sensing protocol where the surface spins act as
reporter spins to sense nuclear spins on the surface.

The nitrogen-vacancy (NV) center in diamond is a promising technology for
quantum sensing at the nanoscale due to its localization to atomic dimensions,
remarkable chemical stability, and unique photophysics Maletinsky  _et al._
(2012); Mamin  _et al._ (2013); Rondin  _et al._ (2014); Müller  _et al._
(2014). The NV center is composed of a substitutional nitrogen atom adjacent
to a lattice vacancy. In order to sense external targets, NV centers must be
brought close to the surface. However, structural and chemical defects at the
surface can give rise to paramagnetic species that cause rapid decoherence of
shallow NV centers Rosskopf  _et al._ (2014); Myers  _et al._ (2014); Romach
_et al._ (2015); Fávaro de Oliveira  _et al._ (2017); Sangtawesin  _et al._
(2019). Some of these defects give rise to unpaired electrons with $S=1/2$ and
a _g_ factor of the free electron Grinolds  _et al._ (2014), and a recent
demonstration utilized this single electron spin as a “quantum reporter” to
detect neighboring proton spins Sushkov  _et al._ (2014). These surface spins
can be stable in ambient conditions over many days Sushkov  _et al._ (2014).
The double electron-electron resonance (DEER) signal associated with surface
spins is diminished after annealing the diamond at 465 in an O2 atmosphere Chu
_et al._ (2014); Lovchinsky  _et al._ (2016); Sangtawesin  _et al._ (2019).
Interestingly, this temperature is believed to remove $sp^{2}$ hybridized
carbon and result in a more perfect oxygen terminated surface Fu  _et al._
(2010); Sangtawesin  _et al._ (2019). Furthermore, the DEER signal is also
suppressed after annealing in vacuum at 650 at which temperature a surface
dangling bond (DB) defect is also removed, as evident in X-ray absorption
Dwyer  _et al._ (2021). These observations suggest that the surface spins
might correspond to $sp^{3}$ DB species. Since $sp^{3}$ DB species can
introduce deep acceptor states into the band gap of diamond they may be a
source of charge and spin fluctuation in near-surface NV centers, which may be
activated even in dark without any illumination Chou  _et al._ (2018);
Bluvstein  _et al._ (2019); Yuan  _et al._ (2020). However, naive models of
dangling bonds at the (100) diamond surface would not be expected to be air-
stable. DEER measurements Sushkov  _et al._ (2014) could not identify the
exact chemical environment of hydrogen atoms at the surface, and we further
note that once Carr-Purcell-Meiboom-Gill protocols are applied in the
detection scheme of the protons, signals associated with 1H and 13C nuclear
spins can be overlapping Loretz  _et al._ (2015), leading to further ambiguity
in extracting direct information about the structure of the diamond surface. A
clear understanding of the diamond surface spin will be crucial for further
improving the sensitivity and resolution of diamond nanosensors.

In this Letter, we propose a simple model for the near-surface DB defect on
(100) oxygenated diamond surface, and we characterize some of their properties
by density functional theory (DFT) calculations. We create a single $sp^{3}$
carbon DB defect on the trench site with a terminating hydroxyl (OH) radical,
see Fig. 1, that can naturally occur at step edges of a polished (100) diamond
surface. The proposed annihilation mechanism of these surface spins is OH
thermal desorption, which induces minimum surface reconstruction. We estimate
the temperature-dependent desorption rate of OH by their DFT activation energy
barriers, and we find correspondence between the annealing of $sp^{3}$ DBs in
the experiments and the simulated OH desorption. Finally, we calculate the
hyperfine structure of four surface spin models, which provides direct
information about the interaction between surface spins and nuclear spins as a
tool for supporting our model and is relevant for reporter spin protocols for
sensing nuclear spins.

First-principles calculations based on DFT were performed using the plane-wave
based Vienna _Ab Initio_ Simulation package Kresse and Furthmüller (1996). The
interactions between the ions and valence electrons were treated by the
projector augmented-wave Blöchl (1994) method. Constant volume relaxation
using a cutoff energy of 370(740) eV in the plane-wave expansion for the wave
function (charge density) resulted in an equilibrium lattice parameter of 3.57
Å, which is only 0.08% larger than the experimental value of 3.567 Å. For the
Brillouin-zone integration of a (2$\times$1) surface unit cell we use a
4$\times$8$\times$1 grid in the Monkhorst-Pack scheme Monkhorst and Pack
(1976). The thickness of vacuum layer is 10 Å and the DBs on the surface are
saturated by hydrogen atoms. We have allowed for all the atoms to fully relax
until the forces were below 0.01 eV/Å except for those at the bottom four
layers which are fixed at their respective position of the (2$\times$1)
reconstructed surface. Afterward we use a (6$\times$6) supercell to simulate a
step diamond model with a surface spin. Here, we use mixed H/O/OH radicals
terminating the diamond surface to minimize the surface states intrusion into
the band gap Kaviani  _et al._ (2014). We apply $\Gamma$-point to map the
Brillouin-zone, and we have tested that it is sufficient to ensure a good
convergence in the total energy difference. The energy barriers of OH
desorption are determined by using nudged elastic band (NEB) method Henkelman
_et al._ (2000). For the exchange-correlation functionals, we used gradient-
corrected Perdew-Burke-Ernzerhof (PBE) functional Perdew  _et al._ (1996) in
the geometry optimization and NEB procedures. Our previous studies of defects
in diamond Gali  _et al._ (2009); Deák  _et al._ (2010, 2014) show that,
unlike PBE functional, the Heyd-Scuseria-Ernzerhof (HSE) hybrid functional
Heyd  _et al._ (2003) provides correct defect levels and defect-related
electronic transitions within $\sim$0.1 eV to experiments, accurate hyperfine
tensors Szász  _et al._ (2013) and energy barriers Deák  _et al._ (2008), so
we apply HSE functional to calculate these properties. The acceptor level of
the defect is obtained from the energy difference between relaxation energy
and vertical ionization energy which avoids the need of charge correction
Stacey  _et al._ (2018). We carried out these calculations at PBE optimized
coordinates and saddle points as the HSE calculations are computationally
challenging in our large slab model.

We have collected the following facts of the surface spins found in
experimental results Sushkov  _et al._ (2014); Chu  _et al._ (2014);
Lovchinsky  _et al._ (2016); Bluvstein  _et al._ (2019); Sangtawesin  _et al._
(2019); Yuan  _et al._ (2020); Dwyer  _et al._ (2021): (i) the surface spin
correlates with X-ray signatures of $sp^{3}$ DB, (ii) it is a spin 1/2
electronic defect, (iii) it must be located within a few layers of the
surface, and (iv) it is chemically stable in air. We pre-examined about eighty
plausible models (see Supplemental Materials in Ref. sup, ) and we present
here the most likely candidate that is consistent with all available
experimental observations.

In experiment, the surface layer is generally damaged and highly strained
because of mechanical polishing Luo  _et al._ (2021), and this strain can be
mitigated by etching and surface processing Sangtawesin  _et al._ (2019).
Systematic simulation and analysis of such disorganized diamond surface is
computationally prohibitive at _ab initio_ level needed for accurate
electronic structure calculation. Instead, we attempted to set up a simple but
still relevant model to identify the surface spins. Our idea is organized as
the following steps: (i) We start with an atomically flat (100) diamond
surface. It is well-known that the dangling bonds at the surface are highly
reactive, thus the surface spins should be located beneath the topmost
surface. (ii) For all carbon atoms of a pristine (100) diamond surface, there
are two bonds upward bonding and two bonds downward bonding. Removing/adding
one carbon atom always generates two DBs. It is therefore not possible to
construct a single $sp^{3}$ DB with a small number of atoms at the (100)
surface. However, at the diamond (111) surface, each carbon atom possesses
three bonds downward(upward) and one bond upward(downward), thus it is
possible to form a single $sp^{3}$ DB defect within minimum change of the
diamond lattice. Therefore, to form a $sp^{3}$ DB beneath (100) surface, a
(111) facet is essential. (iii) As shown in Fig. 1(a), one can find the tilted
(111) facet residing at the trench site of (100) diamond surface caused by a
single atomic step layer. Reflection anisotropy spectroscopy experimental and
theoretical studies Schwitters  _et al._ (2011) show that single layer step on
high-quality atomically smooth H/C(100)-2$\times$1 surfaces predominantly
follows the Chadi step model Chadi (1987). Thus, we modified the atomistic
flat model and introduced a single atomic Chadi-type step on (100) surface to
increase the complexity of the surface at minimum level and simultaneously to
ensure the presence of “disorder” producing a $sp^{3}$ DB. (iv) While one
carbon at trench site (symbol as “*”) is raised up toward (111) direction,
consequently, a single $sp^{3}$ DB defect is formed at the third layer (symbol
as “+”). Meanwhile, a DB at the floating C(*) atom faces the vacuum and could
be passivated by another atom or molecule in the environment. We emphasize
that the polished oxygenated diamond surface should be disordered to some
degree. Our model represents a _locally disordered_ configuration on (100)
diamond surface which has been found by surface spectroscopy Sangtawesin  _et
al._ (2019). (v) As mentioned above, one additional DB generated by the
floating carbon C(*) can be easily saturated by species in atmosphere. The
most probable one is OH radical. X-ray photoelectron spectroscopy and high
resolution electron energy loss spectroscopy measurements Gao  _et al._
(2008); Akhvlediani  _et al._ (2010) indicate that water molecules can
dissociate and adsorb on hydrogenated diamond surfaces at room temperature.
These dissociatively adsorbed water molecules prefer to aggregate on the
trench site Manelli  _et al._ (2010), thus the floating C(*) atom is likely
passivated by H or OH. The concentration of C-OH on water-rich (100) diamond
surface is obviously reduced after annealing to 600 Gao  _et al._ (2008);
Akhvlediani  _et al._ (2010), thus there is good reason to believe that OH
desorption would be the most possible reaction for surface spin elimination at
elevated temperatures. Moreover, after removal of the OH from the surface, no
significant change will happen on diamond surface. Therefore, in the absence
of obvious surface reconstruction, we assume that the floating C(*) atom is
passivated by OH, as Fig. 1(b), and the desorption of OH is the most plausible
mechanism for surface spin elimination reaction, which would not alter
significantly the surface morphology. (vi) We construct a C(100)-6$\times$6
Chadi step model of diamond surface, as shown in Fig. 1, to set up a defect
forming one DB beneath the top surface. In this model, the surface DBs on
terrace are saturated with H/O/OH groups and the DBs on the trench are
saturated with OH group. The surface spin density in this supercell (area of
15.15$\times$15.15 Å2) is 4.4$\times$1013 $\mu$B/cm2 which is consistent with
experiment estimation of $10^{12}\sim 10^{13}\mu$B/cm2 Rosskopf  _et al._
(2014); Myers  _et al._ (2014), although, other defect spins than sp3 DBs may
be also present in (100) oxygenated diamond surface Stacey  _et al._ (2018).
Furthermore, our model implies that even small molecules, e.g., O2 and H2,
cannot penetrate into the surface step and terminate the DB. The $sp^{3}$ DB
apparently can be well protected by the top surface of oxygenated diamond
under ambient conditions.

In our proposed model, the surface spin defect is located three atomic layers
beneath the (100) diamond. Because the DB defect is very close the surface and
the step edge is a chemically reactive site Manelli  _et al._ (2010), the
local surroundings might have the ability of altering the physical properties
or structural stability of the surface DB defect. Therefore we consider three
models, depicting different terminators around the step edge sites. In these
models, all C(*) atoms are saturated with OH and the total spin of these
systems is $S=1/2$. As shown in Fig. 2(a), (1) O/H/H model, one oxygen is
located on the step bridge site and two trench carbon atoms are saturated by
hydrogen atoms; (2) O/OH/OH model, one oxygen is located on the step bridge
site and two trench carbon atoms are saturated by hydroxyl radicals; (3) OH/OH
model, two trench carbon atoms are saturated by hydroxyl radicals. Before
studying the temperature influence for OH desorption, it is important to have
a detailed picture of the electronic properties of these surface spin defect
models.

The HSE results of O/H/H model before and after OH desorption are shown in
Fig. 2(b). The filled(empty) $sp^{3}$ DB levels are located at around the
position of 0.30(3.3) eV above valence band maximum ($E_{\mathrm{VBM}}$).
After OH desorption, DB states are passivated and only some surface states
remain in the bandgap at the position of about 0.7 eV below conduction band
minimum (CBM) which are surface C-H and C-OH image states Kaviani  _et al._
(2014). The calculated energy levels of these three models are shown in Fig.
2(c). The occupied DB states of these three models are at
$E_{\mathrm{VBM}}+0.30$, $+0.22$, $+0.36$ eV, respectively, whereas the empty
level scatters around $E_{\mathrm{VBM}}+3.3\pm 0.07$ eV. As can be seen the
adjacent terminators make only minor changes to the energy level of $sp^{3}$
DB states. We note that the absolute values of these levels are subject of
supercell size effect and they shift down by about 0.17 eV sup . The
calculated vertical ionization energy is around $E_{\mathrm{VBM}}+2.1$ eV but
the ions significantly relax upon ionization leading to around $0.7$ eV
relaxation energy. The final calculated acceptor level of DB defect is around
$E_{\mathrm{VBM}}+1.4$ eV in the $(6\times 6)$ supercell which may have a
larger than $0.1$ eV uncertainty because of the supercell size effect noted
above.

The electron affinity of the oxygenated diamond model is slightly positive as
0.5 eV Kaviani  _et al._ (2014) represents that adsorption of electron
acceptor molecule in atmosphere, e.g., water, will cause upward band bending
Ristein  _et al._ (2004). In realistic oxygenated diamonds this effect may be
significant. Recent experiments could map the band bending in (100) oxygenated
diamond surface Broadway  _et al._ (2018) which indeed found a substantial
upward band bending under ambient conditions. According to their findings, the
quasi Fermi-level position at the diamond surface is around
$E_{\mathrm{VBM}}+0.54$ eV. Taking all the uncertainties in the calculated
acceptor levels of $sp^{3}$ DBs into account one can safely conclude that the
$sp^{3}$ DBs exist in their neutral charge state possessing $S=1/2$ spin
configuration on the diamond surface under realistic conditions.

The $sp^{3}$ DBs can be eliminated by high temperature annealing in O2
atmosphere. It is computationally prohibitive to model the desorption process
in O2 gas, however, it should be possible to carry out annealing in vacuum
that can be directly modeled at _ab initio_ level. We propose a spin
annihilation mechanism that involves OH desorption from the spin defect site.
To verify our argument, we calculate the OH desorption rate and estimate the
number of desorbed OH molecules in a limited area under different
temperatures. To calculate the desorption rate, we introduce the Polanyi-
Wigner equation King (1975),

$R_{\text{des}}=-d\theta/dt=\nu_{n}\exp(-E_{\text{des}}/k_{\text{B}}T)\theta_{n}\text{,}$
(1)

where $E_{\text{des}}$ is the desorption energy barrier, $k_{\text{B}}$ is the
Boltzmann constant, $T$ is temperature, $\nu_{n}$ is the frequency prefactor
which is typically 1015 1/s due to the larger entropy of the desorbed
molecules in vacuum Alfé and Gillan (2006, 2007). $\theta$ is the coverage of
OH on DB site, which is equal to one here. Because the activation energy
barriers of OH desorption would possibly be influenced by the local
environment, we calculate the energy barriers of the three considered models
as Fig. 2(a). An illustration of OH thermal desorption from diamond step edge
is shown in Fig. 3(a). The calculated energy barriers in these three models
are 0.89, 1.12, and 0.96 eV, respectively. Next, we use Eq. 1 to calculate the
OH desorption rates as a function of temperature, as shown in Fig. 3(b). The
results indicate that the number of desorbed OH groups (in one hour) at 600 is
at least one order of magnitude more than that at 465 . Compared with the
experimental estimation of surface density ($10^{12}\sim 10^{13}$ $\mu$B/cm2),
the desorption energies in the range of 0.89 $\sim$ 1.12 eV result in the
removal of such amount of surface spins on a surface area of 1 cm2 at 600 in
one hour. At lower temperature of 465 it requires more than three hours to
remove all spins in the same area which is consistent with the experiment
observations Chu  _et al._ (2014); Lovchinsky  _et al._ (2016).

So far we have addressed the microscopic structure of diamond surface spins.
With a single electron occupying the near-surface DB the net electron spin is
$S=1/2$, which means this near-surface DB spin can interact with the present
13C in diamond or 1H isotopes on the surface (nuclear spin $I=1/2$). The
hyperfine interaction (HFI) between electron spin S and a set of nuclear spins
Il ($l$=1,2…,N) can be described by the Hamiltonian H=$\sum$ SA(l)Il with A(l)
being the hyperfine coupling tensor associated with the $l$th nucleus at site
Rl. In atomic units the hyperfine tensor can be written as:
$A_{ij}^{(l)}=\frac{8\pi}{3}g_{e}\mu_{e}g_{n}\mu_{n}\rho_{n}(\textbf{{R}}_{l})+g_{e}\mu_{e}g_{n}\mu_{n}\int
d\textbf{{r}}\frac{3r_{i}r_{j}-r^{2}\delta_{ij}}{r^{5}}\rho_{s}(\textbf{{{r}}})$.
The first term ($a_{iso}$) provides the isotropic hyperfine interaction and is
referred to as the Fermi-contact term, where $g_{e}(g_{n})$ is the
gyromagnetic ratio of electron(nucleus), $\mu_{e}$($\mu_{n}$) is the Bohr
magneton of electron(nucleus), and $\rho_{n}$($\rho_{s}$) is the
nuclear(electron) spin density. This term is proportional to the magnitude of
electron spin density at the nucleus center. The second term provides the
anisotropic HFI and is referred to as the dipolar term where r is the vector
connecting the electron and the nuclear spin. The electron-nucleus Hamiltonian
can be put in a simpler form by letting the magnetic field and the crystal _C_
axis of symmetry lie in the xz plane, where the C axis is in the z direction.
The transformed Hamiltonian is then simplified to
$\textbf{{H}}=a\textbf{{S}}_{Z}\textbf{{I}}_{Z}+b\textbf{{S}}_{Z}\textbf{{I}}_{X}$
Rowan  _et al._ (1965); Schweiger and Jeschke (2001) , where the z axis is
along the applied magnetic field. The quantities $a=A_{ZZ}$ and
$b=(A^{2}_{ZX}+A^{2}_{ZY})^{1/2}$ describe the secular and pseudo-secular
hyperfine couplings. Non-secular terms with $\textbf{{S}}_{X}$ and
$\textbf{{S}}_{Y}$ are neglected. Further, the separation (r) between the
electron spin and nuclear spin and the angle ($\theta$) that the vector
between them makes with applied magnetic field can be extracted from the
fitting of HFI parameters ($a$, $b$), where
$a=a_{iso}+T(3\cos^{2}\theta$$-$$1)$, $b=3T\sin{\theta}\cos{\theta}$, and
$T=g_{e}\mu_{e}g_{n}\mu_{n}/r^{3}$ Schweiger and Jeschke (2001). From the
above description, it is possible to detect the relative position and angle of
a proton around the surface spin.

Combined with initialization and readout at a proximal NV center, our detailed
model for the hyperfine coupling between the dangling bond and proximate
nuclear spins can be used to probe our model for the dangling bond and to do
spectroscopy on surface species. As an example, we start with the extracted
relative positions and angles of two NVs with dark electron spins in Ref.
Sushkov  _et al._ , 2014, labeled as NVA and NVB. We perform precise
calculations of the HFI parameters ($a$,$b$) for the surface spin model. By
scanning the calculated HFI parameters through all C and H atoms, more than
ten atoms possess large values (> 10 MHz), but mostly the $b$ values are small
(< 5 MHz), in contrast to the usual assumptions that the Fermi-contact term is
negligible. As shown in Fig. 2(d), the spin density is mostly localized on the
DB site where it reflects a large Fermi-contact value of (329, 105) MHz on DB
13C atom and a small value of (4.3, 2.2) MHz on surface hydrogen atom of the
OH group. The latter is in good agreement with the extracted $a$ and $b$ for
NVB, (4.0, 2.2) MHz. However, we are unable to find a surface model to explain
the measured HFI parameters for NVA with hyperfine couplings of around 10 MHz
associated with the proton spins. Further investigations using DEER
spectroscopy can elucidate alternate possible surface spin configurations,
surface proton configurations, and subsurface dark spins. Our results
demonstrate that the combination of DEER spectrum and _ab initio_ simulations
is necessary in the reporter spin protocol for sensing of nuclear spins
because the Fermi-contact term can be sizable and no simple approximation is
able to estimate its strength.

In conclusion, by means of DFT calculations, we propose a reliable model of
surface spins on the (100) diamond surface. A (111) facet is essential to
create a DB defect beneath diamond surface and this facet naturally exists at
the step of (100) diamond surface. We also demonstrated that the OH desorption
annihilates the surface spin. HFI calculations could provide a great help in
the identification of the defect structure in detail. Our present work
indicates large HFI 13C parameters of ($a$, $b$)=(337, 106) MHz for the
surface dangling bond. The identification and characterization of surface
spins on diamond provided a strong hint to continue engineering of high-
quality materials for device applications. A new avenue is to extend these
studies to (111) diamond surface where step-free surface can be grown Tokuda
_et al._ (2012) with proper surface termination.

We acknowledge that the results of this research have been achieved using the
DECI resource Eagle HPC based in Poland at Poznan with support from the PRACE
aisbl and resources provided by the Hungarian Governmental Information
Technology Development Agency. A. G. acknowledges the National Research,
Development, and Innovation Office of Hungary grant No. KKP129866 of the
National Excellence Program of Quantum-coherent materials project and the
Quantum Information National Laboratory supported by the Ministry of
Innovation and Technology of Hungary, as well as the EU Commission for the
H2020 Quantum Technology Flagship project ASTERIQS (Grant No. 820394). NPdL
was supported by the DARPA DRINQS program (grant D18AC00015). J. C.
acknowledges the financial support from the Ministry of Science and
Technology, Taiwan (MOST 109-2112-M-018-008-MY3).

## References

  * Maletinsky  _et al._ (2012) P. Maletinsky, S. Hong, M. S. Grinolds, B. Hausmann, M. D. Lukin, R. L. Walsworth, M. Loncar, and A. Yacoby, Nature Nanotechnology 7, 320 (2012).
  * Mamin  _et al._ (2013) H. J. Mamin, M. Kim, M. H. Sherwood, C. T. Rettner, K. Ohno, D. D. Awschalom, and D. Rugar, Science 339, 557 (2013).
  * Rondin  _et al._ (2014) L. Rondin, J.-P. Tetienne, T. Hingant, J.-F. Roch, P. Maletinsky, and V. Jacques, Reports on Progress in Physics 77, 056503 (2014).
  * Müller  _et al._ (2014) C. Müller, X. Kong, J.-M. Cai, K. Melentijević, 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, Nature Communications 5, 4703 (2014).
  * Rosskopf  _et al._ (2014) T. Rosskopf, A. Dussaux, K. Ohashi, M. Loretz, R. Schirhagl, H. Watanabe, S. Shikata, K. Itoh, and C. Degen, Physical Review Letters 112, 147602 (2014).
  * Myers  _et al._ (2014) B. Myers, A. Das, M. Dartiailh, K. Ohno, D. Awschalom, and A. Bleszynski Jayich, Physical Review Letters 113, 027602 (2014).
  * Romach  _et al._ (2015) 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, Physical Review Letters 114, 017601 (2015).
  * Fávaro de Oliveira  _et al._ (2017) F. Fávaro de Oliveira, D. Antonov, Y. Wang, P. Neumann, S. A. Momenzadeh, T. Häßermann, A. Pasquarelli, A. Denisenko, and J. Wrachtrup, Nature Communications 8, 15409 (2017).
  * Sangtawesin  _et al._ (2019) S. Sangtawesin, B. L. Dwyer, S. Srinivasan, J. J. Allred, L. V. H. Rodgers, K. De Greve, A. Stacey, N. Dontschuk, K. M. O’Donnell, D. Hu, D. A. Evans, C. Jaye, D. A. Fischer, M. L. Markham, D. J. Twitchen, H. Park, M. D. Lukin, and N. P. de Leon, Phys. Rev. X 9, 031052 (2019).
  * Grinolds  _et al._ (2014) M. S. Grinolds, M. Warner, K. D. Greve, Y. Dovzhenko, L. Thiel, R. L. Walsworth, S. Hong, P. Maletinsky, and A. Yacoby, Nature Nanotechnology 9, 279 (2014).
  * Sushkov  _et al._ (2014) A. Sushkov, I. Lovchinsky, N. Chisholm, R. Walsworth, H. Park, and M. Lukin, Physical Review Letters 113, 197601 (2014).
  * Chu  _et al._ (2014) Y. Chu, N. de Leon, B. Shields, B. Hausmann, R. Evans, E. Togan, M. J. Burek, M. Markham, A. Stacey, A. Zibrov, A. Yacoby, D. Twitchen, M. Loncar, H. Park, P. Maletinsky, and M. Lukin, Nano Letters 14, 1982 (2014).
  * Lovchinsky  _et al._ (2016) I. Lovchinsky, A. O. Sushkov, E. Urbach, N. P. d. Leon, S. Choi, K. D. Greve, R. Evans, R. Gertner, E. Bersin, C. Müller, L. McGuinness, F. Jelezko, R. L. Walsworth, H. Park, and M. D. Lukin, Science 351, 836 (2016).
  * Fu  _et al._ (2010) K.-M. C. Fu, C. Santori, P. E. Barclay, and R. G. Beausoleil, Applied Physics Letters 96, 121907 (2010).
  * Dwyer  _et al._ (2021) B. L. Dwyer, L. V. H. Rodgers, E. K. Urbach, D. Bluvstein, S. Sangtawesin, H. Zhou, Y. Nassab, M. Fitzpatrick, Z. Yuan, K. De Greve, E. L. Peterson, J.-P. Chou, A. Gali, V. V. Dobrovitski, M. D. Lukin, and N. P. de Leon, arXiv:2103.12757 [cond-mat, physics:quant-ph] (2021), arXiv: 2103.12757.
  * Chou  _et al._ (2018) J.-P. Chou, Z. Bodrog, and A. Gali, Physical Review Letters 120 (2018), 10.1103/PhysRevLett.120.136401.
  * Bluvstein  _et al._ (2019) D. Bluvstein, Z. Zhang, and A. C. B. Jayich, Physical Review Letters 122 (2019), 10.1103/PhysRevLett.122.076101.
  * Yuan  _et al._ (2020) Z. Yuan, M. Fitzpatrick, L. V. H. Rodgers, S. Sangtawesin, S. Srinivasan, and N. P. de Leon, Physical Review Research 2, 033263 (2020), publisher: American Physical Society.
  * Loretz  _et al._ (2015) M. Loretz, J. Boss, T. Rosskopf, H. Mamin, D. Rugar, and C. Degen, Physical Review X 5, 021009 (2015).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Physical Review B 54, 11169 (1996).
  * Blöchl (1994) P. E. Blöchl, Physical Review B 50, 17953 (1994).
  * Monkhorst and Pack (1976) H. J. Monkhorst and J. D. Pack, Physical Review B 13, 5188 (1976).
  * Kaviani  _et al._ (2014) M. Kaviani, P. Deák, B. Aradi, T. Frauenheim, J.-P. Chou, and A. Gali, Nano Letters 14, 4772 (2014).
  * Henkelman  _et al._ (2000) G. Henkelman, B. P. Uberuaga, and H. Jònsson, The Journal of Chemical Physics 113, 9901 (2000).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Physical Review Letters 77, 3865 (1996).
  * Gali  _et al._ (2009) A. Gali, E. Janzén, P. Deák, G. Kresse, and E. Kaxiras, Physical Review Letters 103, 186404 (2009).
  * Deák  _et al._ (2010) P. Deák, B. Aradi, T. Frauenheim, E. Janzén, and A. Gali, Phys. Rev. B 81, 153203 (2010).
  * Deák  _et al._ (2014) P. Deák, B. Aradi, M. Kaviani, T. Frauenheim, and A. Gali, Physical Review B 89, 075203 (2014).
  * Heyd  _et al._ (2003) J. Heyd, G. E. Scuseria, and M. Ernzerhof, The Journal of Chemical Physics 118, 8207 (2003), https://doi.org/10.1063/1.1564060 .
  * Szász  _et al._ (2013) K. Szász, T. Hornos, M. Marsman, and A. Gali, Physical Review B 88, 075202 (2013).
  * Deák  _et al._ (2008) P. Deák, B. Aradi, T. Frauenheim, and A. Gali, Materials Science and Engineering: B 154-155, 187 (2008).
  * Stacey  _et al._ (2018) A. Stacey, N. Dontschuk, J.-P. Chou, D. A. Broadway, A. K. Schenk, M. J. Sear, J.-P. Tetienne, A. Hoffman, S. Prawer, C. I. Pakes, A. Tadich, N. P. de Leon, A. Gali, and L. C. L. Hollenberg, Advanced Materials Interfaces 6, 1801449 (2018).
  * (33) See Suplemental Material.
  * Luo  _et al._ (2021) H. Luo, K. M. Ajmal, W. Liu, K. Yamamura, and H. Deng,  3, 022003 (2021), publisher: IOP Publishing.
  * Schwitters  _et al._ (2011) M. Schwitters, D. S. Martin, P. Unsworth, T. Farrell, J. E. Butler, M. Marsili, O. Pulci, and P. Weightman, Physical Review B 83, 085402 (2011).
  * Chadi (1987) D. J. Chadi, Physical Review Letters 59, 1691 (1987).
  * Gao  _et al._ (2008) X. Gao, L. Liu, D. Qi, S. Chen, A. T. S. Wee, T. Ouyang, K. P. Loh, X. Yu, and H. O. Moser, The Journal of Physical Chemistry C 112, 2487 (2008).
  * Akhvlediani  _et al._ (2010) R. Akhvlediani, S. Michaelson, and A. Hoffman, Surface Science 604, 2129 (2010).
  * Manelli  _et al._ (2010) O. Manelli, S. Corni, and M. C. Righi, The Journal of Physical Chemistry C 114, 7045 (2010).
  * Ristein  _et al._ (2004) J. Ristein, M. Riedel, and L. Ley, Journal of The Electrochemical Society 151, E315 (2004).
  * Broadway  _et al._ (2018) D. A. Broadway, N. Dontschuk, A. Tsai, S. E. Lillie, C. T.-K. Lew, J. C. McCallum, B. C. Johnson, M. W. Doherty, A. Stacey, L. C. L. Hollenberg, and J.-P. Tetienne, Nature Electronics 1, 502 (2018).
  * King (1975) D. A. King, Surface Science 47, 384 (1975).
  * Alfé and Gillan (2006) D. Alfé and M. J. Gillan, Journal of Physics: Condensed Matter 18, L451 (2006).
  * Alfé and Gillan (2007) D. Alfé and M. J. Gillan, The Journal of Chemical Physics 127, 114709 (2007).
  * Rowan  _et al._ (1965) L. G. Rowan, E. L. Hahn, and W. B. Mims, Physical Review 137, A61 (1965).
  * Schweiger and Jeschke (2001) A. Schweiger and G. Jeschke, (Oxford University Press, 2001).
  * Tokuda  _et al._ (2012) N. Tokuda, T. Makino, T. Inokuma, and S. Yamasaki, Japanese Journal of Applied Physics 51, 090107 (2012).

Figure 1:  (Color online) (a) Top and side view of step model of
C(100)-6$\times$6 diamond surface. The (100) and tilted (111) facets are
highlighted with red and blue colors, respectively. (b) Top view of surface
spin model. The DB position is labeled with symbol “+”. Figure 2:  (Color
online) (a) Top view of three different surface spin defect models. The color
settings are the same as those in Fig. 1. (b) The HSE results of the energy
levels of O/H/H model before and after OH desorption. (c) The HSE energy
levels of three surface spin defect models. The valence and conduction bands
are depicted as blue and red regions, respectively. The valence band maximum
is aligned to zero. (d) The isosurface of the calculated spin density
(isovalues are 5$\times$102 e/Bohr3 in main figure and 1.3$\times$104 e/Bohr3
in the top-left inset) for surface spin model. The spin density is clearly
localized on the DB site. The carbon, oxygen, and hydrogen atoms are grey,
red, and pink balls, respectively. Figure 3:  (Color online) (a)The schematic
illustration of surface spin (green lobes) annihilation. After thermal
annealing in vacuum, the OH radical desorbs from the surface trench site,
causing surface reconstruction, and DB is eventually passivated. The color
settings are the same as those in Fig. 1. (b) Calculated OH desorption rate of
three models as a function of temperature. The annihilation temperatures of
465 and 600 are denoted as dashed and dash-dotted lines, respectively.