# Online Active Learning for Soft Sensor Development using Semi-Supervised
Autoencoders

Davide Cacciarelli<EMAIL_ADDRESS>
Murat Kulahci<EMAIL_ADDRESS>
Department of Applied Mathematics and Computer Science
Technical University of Denmark
Kgs. Lyngby, Denmark  and John Tyssedal<EMAIL_ADDRESS>
Department of Mathematical Sciences
Norwegian University of Science and Technology
Trondheim, Norway

###### Abstract

Data-driven soft sensors are extensively used in industrial and chemical
processes to predict hard-to-measure process variables whose real value is
difficult to track during routine operations. The regression models used by
these sensors often require a large number of labeled examples, yet obtaining
the label information can be very expensive given the high time and cost
required by quality inspections. In this context, active learning methods can
be highly beneficial as they can suggest the most informative labels to query.
However, most of the active learning strategies proposed for regression focus
on the offline setting. In this work, we adapt some of these approaches to the
stream-based scenario and show how they can be used to select the most
informative data points. We also demonstrate how to use a semi-supervised
architecture based on orthogonal autoencoders to learn salient features in a
lower dimensional space. The Tennessee Eastman Process is used to compare the
predictive performance of the proposed approaches.

Keywords: Active Learning, Semi-Supervised Learning, Linear Regression,
Autoencoder.

##  1 Introduction

In industrial operations, soft sensors are frequently used for real-time
prediction of hard-to-measure process variables, as well as to support system
backup strategies, what-if analysis, sensor validation, and fault diagnosis
(Fortuna et al., 2007). Soft sensors are classified into two types: model-
driven sensors, which are used to depict the ideal steady-state of a process
under normal operating conditions, and data-driven sensors, which are used to
better approximate real process conditions (Kadlec et al., 2009). Many labeled
observations are required for training the regression models used in soft
sensor development, but in industrial contexts, data is often abundant only in
an unlabeled form. Obtaining product information can be both expensive and
time consuming, as it may necessitate the intervention of a human expert or
the use of expensive testing equipment. As a result, active learning is
becoming increasingly useful for reducing the number of labels required to
achieve compelling predictive performance. Active learning-based sampling
schemes use some evaluation criteria to assess the informativeness of the
unlabeled data points and prioritize the labeling of the most useful instances
for building the model. Three macro scenarios can be identified depending on
how the unlabeled instances are fed into the learner and then selected to be
labeled by an oracle (Settles, 2009). The first scenario is referred to as
membership query synthesis, and it allows the learner to query the labels of
synthetically generated instances rather than those sampled from the process
distribution. The second scenario is stream-based active learning, also known
as selective sampling. It denotes a situation in which instances are drawn
sequentially and the learner must immediately decide whether to keep the
instance and query its label or discard it. The third and final scenario is
pool-based active learning, which depicts a situation where a large amount of
unlabeled data is collected all at once and made available to the learner,
which can rank all of the data points and select the most informative ones.
While many researchers have been working on active learning in the latest
years, pool-based active learning for classification has received the most
attention (Cai et al., 2013).

In this work, we focus on stream-based active learning (Cacciarelli and
Kulahci, 2023), which represents a more difficult task as the learner cannot
observe all of the available observations before deciding which labels to
query. We believe that this scenario accurately reflects high-volume
production processes in which samples are processed very rapidly and labels
are no longer retrievable. Stream-based active learning should be considered
and prioritized for all industrial processes with similar properties.

##  2 Background

In regression modeling, we try to learn a function
$\hat{f}:\boldsymbol{x}\in\mathbb{R}^{p}\rightarrow y\in\mathbb{R}$ to predict
a quality characteristic or a hard-to-measure variable $y\in\mathbb{R}$ that
is related to other process variables $\boldsymbol{x}\in\mathbb{R}^{p}$.
Accordingly with many active learning approaches (Cai et al., 2013), we assume
a labeled dataset
$\mathcal{L}=\\{(\boldsymbol{x}_{i},y_{i})_{i=1}^{n},\boldsymbol{x}_{i}\in\mathbb{R}^{p},y\in\mathbb{R}\\}$
with $n$ observations is initially available to fit a linear regression model
of the kind

$f(\boldsymbol{x};\boldsymbol{\beta})=\sum_{i=0}^{p}\beta_{i}x_{i}=\boldsymbol{\beta}^{T}\boldsymbol{x}$

where $x_{0}=1$ is the intercept term and $x_{i}$ with $i=1,...,p$ are the $p$
process variables. Parameters are estimated by minimizing a squared error loss
given by

$\hat{\epsilon}=\frac{1}{n}\sum_{i=1}^{n}(y_{i}-f(\boldsymbol{x}_{i}))^{2}$
(1)

After an initial model has been built, we aim to acquire additional
observations by evaluating the unlabeled data points, until a budget
constraint is met. Some commonly encountered approaches are presented below.

###  2.1 Mahalanobis Distance

If we only examine the feature space, a desirable property that we might
pursue when collecting instances for our training set $\mathcal{L}$, is to
ensure diversity among the observations. The Hotelling $T^{2}$ control chart,
which is widely used in statistical process control (SPC) to detect anomalous
data points (Hotelling, 1947), can accomplish so. As we do in SPC, in this
case we use Mahalanobis distance (Equation 2) to measure the dissimilarity
between the new unlabeled instances and the observations in the current
training set $\mathcal{L}$. The Hotelling $T^{2}$ statistic for a new
unlabeled instance $\boldsymbol{x}$ is computed as

$T^{2}(\boldsymbol{x})=(\boldsymbol{x}-\bar{\boldsymbol{x}})^{T}\boldsymbol{S}^{-1}(\boldsymbol{x}-\bar{\boldsymbol{x}})$
(2)

where $\bar{\boldsymbol{x}}$ and $\boldsymbol{S}$ correspond to the sample
mean vector and sample covariance matrix of $\mathcal{L}$, respectively. This
approach has been extended to a principal component regression (PCR) model by
Ge et al. (2014), who proposed a sampling index dependent on the Hotelling
$T^{2}$ statistic and the squared prediction error. In this case, the sampling
function is simply represented by
$argmax_{\boldsymbol{x}}T^{2}(\boldsymbol{x})$.

###  2.2 Query By Committee

While the previous approach only considers the feature space, query by
committee (QBC) tries to evaluate the uncertainty about the response. This
strategy, initially introduced for classification problems, was extended to
regression tasks by Burbidge et al. (1997). The main intuition is that by
building an ensemble of regression models trained on bootstrap replica of the
original training set $B(K)=\\{f_{1},f_{2},...,f_{K}\\}$ we can approximate
the distribution of the predictive variance. Once the ensemble has been built,
we can measure the variance of the predictions made by the committee members
for each unlabeled observation $\boldsymbol{x}$. This variance, or ambiguity,
is computed as

$a(\boldsymbol{x})=\frac{1}{K}\sum_{i=1}^{K}(f_{i}(\boldsymbol{x})-y_{K}(\boldsymbol{x}))^{2}$
(3)

where $y_{K}(\boldsymbol{x})$ is the mean of the predictions made by the
ensemble members. The sampling function then simply becomes
$argmax_{\boldsymbol{x}}a(\boldsymbol{x})$. The key intuition is that if many
models disagree on the label associated with an instance, that instance is an
ambiguous one.

###  2.3 Expected Model Change

Introduced by Cai et al. (2013), expected model change (EMC) suggests querying
the unlabeled example that would cause the maximum change in the current model
parameters, if we knew its label. The model change is measured as the
difference between the current model parameters and the parameters obtained
after fitting the model on the enlarged training set
$\mathcal{L}^{+}=\mathcal{L}\cup(\boldsymbol{x}^{+},y^{+})$. The gradient of
the loss is used to estimate the model change. Considering the augmented
training set $\mathcal{L}^{+}$, the loss function shown in Equation 1 becomes

$\hat{\epsilon}=\frac{1}{n}\sum_{i=1}^{n}(y_{i}-f(\boldsymbol{x}_{i}))^{2}+(y^{+}-f(\boldsymbol{x}^{+}))^{2}$

where the last term, which is hereinafter referred to as
$\ell_{\boldsymbol{x}^{+}}(\boldsymbol{\beta})$, represents the difference
between the loss measured with the model trained on $\mathcal{L}$ and the one
trained on $\mathcal{L}^{+}$. The derivative of the marginal loss
$\ell_{\boldsymbol{x}^{+}}(\boldsymbol{\beta})$ with respect to the parameters
$\boldsymbol{\beta}$ in the new point $\boldsymbol{x}^{+}$ is given by

$\displaystyle\hat{\epsilon}=\frac{\partial\ell_{\boldsymbol{x}^{+}}(\boldsymbol{\beta})}{\partial\boldsymbol{\beta}}$
$\displaystyle=2(y^{+}-f(\boldsymbol{x}^{+}))\frac{\partial
f(\boldsymbol{x}^{+})}{\partial\boldsymbol{\beta}}$ (4)
$\displaystyle=2(y^{+}-f(\boldsymbol{x}^{+}))\frac{\partial\boldsymbol{\beta}^{T}\boldsymbol{x}^{+}}{\partial\boldsymbol{\beta}}$
$\displaystyle=2(y^{+}-f(\boldsymbol{x}^{+}))\boldsymbol{x}^{+}$

Since we do not know the true label of $\boldsymbol{x}^{+}$, $y^{+}$ it is
going to be replaced by the predictions $f_{i}(\boldsymbol{x}^{+})$ made by
the members of the bootstrap ensemble $B(K)$. Finally, the sampling function
is given by
$argmax_{\boldsymbol{x}}\frac{1}{K}\sum_{i=1}^{K}\|(f_{i}(\boldsymbol{x})-f(\boldsymbol{x}))\boldsymbol{x}\|$.

##  3 Proposed Approach

Given the impossibility of ranking unlabeled instances in real-time and
deterministically optimizing the sampling criteria, we propose leveraging
unlabeled data to impose a threshold, or upper control limit (UCL), on the
informativeness of the incoming data points. The unlabeled data pool can be
acquired by either observing the process for a period of time without sampling
the product information $y$ or by using data that is already available in the
form of a historical database
$\mathcal{H}=\\{(\boldsymbol{x}_{i}),\boldsymbol{x}_{i}\in\mathbb{R}^{p}\\}$.
The primary difference between pool-based active learning and online active
learning is that the labels of observations pertaining to $\mathcal{H}$ can no
longer be queried because they only exist digitally, and the associated
physical part or component is no longer available. The data in $\mathcal{H}$
is used to estimate the distribution of the statistics employed by the
criteria in Equations 2, 3, and 4. In this study, we employed kernel density
estimation with a Gaussian kernel. The UCL is then determined by specifying
the appropriate sampling rate $\alpha$. For a given criterion $\mathcal{J}$,
the threshold is defined as

$P(\mathcal{J(\boldsymbol{x})}\geq UCL)=\alpha$ (5)

Using the UCL obtained from Equation 5, we should then only collect the
$\alpha$-percent most informative data points, according to the specific
criterion $\mathcal{J}$ (Cacciarelli et al., 2022, 2023). In this work, we
test the stream-based active learning routine using Mahalanobis distance,
ambiguity, and expected model change as sampling criteria $\mathcal{J}$.
Before starting the active learning routine and collecting additional
observations, we also propose to use a semi-supervised architecture by
training an autoencoder network on the historical data $\mathcal{H}$. With
semi-supervised learning, we can exploit all the available unlabeled data and
learn how to extract relevant features that could be better predictors than
the raw input features. Indeed, when variables are highly correlated, it has
been demonstrated in the literature that a PCR model can be enhanced with
semi-supervision (Frumosu and Kulahci, 2018). With regards to deep learning
methods, autoencoders have been proposed to deal with semi-supervised learning
in fault classification (Jia et al., 2020; Jiang et al., 2017). Recently,
autoencoders have also been investigated in soft sensors applications (Yuan et
al., 2018; de Almeida Moreira et al., 2021) but their contribution to the
stream-based active learning scenario has not been evaluated yet. In this
work, we propose the use of a semi-supervised architecture as the one shown in
Figure 1. An orthogonal autoencoder (OAE) is employed for feature extraction
and the encoded features are then used as predictors in a linear regression
model. An OAE is an autoencoder network that minimizes an Ortho-Loss, which is
comprised of a squared reconstruction loss and an orthogonality regularization
term (Wang et al., 2019; Cacciarelli and Kulahci, 2022). The regularization,
weighted by a parameter $\lambda$, encourages the network to learn
uncorrelated features in its bottleneck. This is particularly beneficial to
alleviate the multicollinearity issue in the regression modeling stage.

Figure 1: Semi-supervised architecture based on OAE.

The main advantage of the semi-supervised model is that the extracted features
are more expressive than the original process variables. However, if the
dimensionality of the bottleneck is lower than the one of the input features,
there is an additional benefit for active learning. Indeed, because the
majority of the provided active learning approaches are model-based, an
initial number of labels is required. QBC and EMC, in particular, employ the
linear regression model’s predictions to select the data points that should be
queried. To uniquely determine the coefficients of a regression model, we need
a number of observations larger than the number of parameters
$\boldsymbol{\beta}$ to be estimated. This initial set of observations is
usually collected at random (Cai et al., 2013). As a result, by reducing the
dimensionality of the parameters $\boldsymbol{\beta}$, we will be able to
anticipate the active collection phase and get more robust estimates for the
same experimental cost. The complete active learning routine is reported
Appendix A.

##  4 Experiments

The Tennessee Eastman Process (TEP) is considered the gold-standard benchmark
for testing process control approaches (Ricker, 1996; Capaci et al., 2019)
and, recently, it has also been used for validating active learning and soft
sensor development methods (Zhu et al., 2015; Grbić et al., 2013). The
variables that have been used as predictors in regression modeling are the
same 16 controlled process variables used by Zhu et al. (2015) and Grbić et
al. (2013), and the continuous response is Stream 9E, a composition
measurement belonging to the purge stream.

(a) Semi-Supervised Learning

(b) Active Learning

Figure 2: (a) shows the learning curves with random sampling using the
original process variables and the features extracted by the OAE and (b) shows
the learning curves of different active learning methods: random (RND),
Hotelling $T^{2}$ (HOT), query by committee (QBC), and expected model change
(EMC).

Figure 2a shows how the suggested semi-supervised architecture can increase
the predictive performance. Data is randomly sampled in both situations, but
the two linear models are fitted using the original process variables and the
features extracted by the OAE. To ensure comparability between the two
learning curves, we fitted the first model when the number of gathered
observations exceeded 16. We believe the improvement is due to the OAE’s
ability to express nonlinear relationships in data in its encoded features and
to the fact that with the extracted features, we have the same number of
observations to estimate a smaller number of parameters. In Figure 2b, we try
to improve the semi-supervised result by using the proposed active learning
strategies. It is clear that EMC and QBC consistently outperform the passive
random approach in terms of recommending the most informative data points. On
the contrary, the Mahalanobis distance appears to worsen the predictive
performance. We believe this may due to the fact that data points with high
$T^{2}$ statistics may be outliers, whose inclusion in the training set
eventually degrades performance. It should be noted that in Figure 2b all the
sampling strategies use the features extracted by the OAE. Experimental setup
and training details are reported in Appendix B.

##  5 Conclusion

Industrial data is often only available unlabeled as quality inspections and
manual annotation tasks are costly and time-consuming. In this work, we
proposed a semi-supervised model based on OAEs for extracting relevant
features and reducing multicollinearity. On top of this, we reviewed and
adapted for the online setting some of the most widely used active learning
strategies for linear regression. The analysis demonstrates how properly using
the historical data and taking into account the expected response allows for a
faster reduction of the prediction error. For future research, we will
consider more advanced architectures such as LSTM autoencoders or transformers
to obtain encoded features that take into account the temporal dependency in
the data.

## Appendix A. Online Active Learning Routine

Algorithm 1

1:a historical unlabeled dataset $\mathcal{H}$, a labeled dataset
$\mathcal{L}$, a data stream $\mathcal{S}$, a budget $b$, and a criterion
$\mathcal{J}$.

2:Train an OAE on $\mathcal{H}$

3:Encode observations in $\mathcal{H}$ and $\mathcal{L}$:
$\boldsymbol{x}\in\mathbb{R}^{p}\longrightarrow\boldsymbol{z}\in\mathbb{R}^{k}$

4:Fit a linear regression model on the encoded features $\boldsymbol{z}$ and
labels $y$ obtained from $\mathcal{L}$

5:Compute $\mathcal{J}$ on the encoded features $\boldsymbol{z}$ pertaining to
$\mathcal{H}$ and estimate a threshold (UCL) using Equation 5

6:$i\leftarrow 0$, $c\leftarrow 0$

7:while $c\leq b$ and $i\leq|\mathcal{S}|$ do

8: Encode ith observation from the stream $\mathcal{S}$:
$\boldsymbol{x}_{i}\in\mathbb{R}^{p}\longrightarrow\boldsymbol{z}_{i}\in\mathbb{R}^{k}$

9: if $\mathcal{J}(\boldsymbol{z})\geq UCL$ then

10: Ask for the label $y_{i}$ and augment the labeled dataset:
$\mathcal{L}^{+}=\mathcal{L}\cup(\boldsymbol{z}_{i},y_{i})$

11: $c\leftarrow c+1$

12: Update model (repeat Step 3)

13: Update threshold (repeat Step 4)

14: else

15: Discard $\boldsymbol{x}_{i}$

16: end if

17: $i\leftarrow i+1$

18:end while

## Appendix B. Experimental Setup and Training Details

The data is generated using the MATLAB code provided by Reinartz et al. (2021)
and Andersen et al. (2022) with the Ricker closed-loop simulation model. No
faults have been introduced throughout the 50 simulation runs, which are
generated providing different seeds to the simulator. The variables used are
reported in table 1. Sample rate was set to 1 minute.

The active learning routine was tested once on each of the 50 simulation runs.
Figure 2 reports the mean and standard deviation for each method across these
50 runs (shaded regions indicate $\pm$ 1 standard deviation).

With regards to the autoencoder structure, we used an encoder whose
dimensionality of the layers corresponds to [16, 160, 80, 40, 20, 10]. The
decoder is symmetrical to the encoder. The penalty term corresponding to the
weight of the orthogonality regularization in the loss function was set to
0.10. No fixed number of epochs was used for the training as we followed an
early stopping approach, setting a patience of 10 on the number of accepted
epochs without improvement on the validation loss (20% of the training data is
used for validation). Finally, the bandwidth used for the kernel density
estimation of the UCL is found using Scott’s rule (Scott, 1992).

Process Variable | ID
---|---
A Feed (Stream 1) | XMEAS 1
D Feed (Stream 2) | XMEAS 2
E Feed (Stream 3) | XMEAS 3
A and C Feed (Stream 4) | XMEAS 4
Recycle Flow (Stream 8) | XMEAS 5
Reactor Feed Rate (Stream 6) | XMEAS 6
Reactor Temperature | XMEAS 9
Purge Rate (Stream 9) | XMEAS 10
Separator Temperature | XMEAS 11
Separator Pressure | XMEAS 13
Product Separator Underflow (Stream 10) | XMEAS 14
Stripper Pressure | XMEAS 16
Stripper Temperature | XMEAS 18
Stripper Steam Flow | XMEAS 19
Reactor Cooling Water Outlet Temperature | XMEAS 21
Separator Cooling Water Outlet Temperature | XMEAS 22
Table 1: Monitored variables of the TEP.

## References

  * Andersen et al. (2022) Emil B. Andersen, Isuru A. Udugama, Krist V. Gernaey, Abdul R. Khan, Christoph Bayer, and Murat Kulahci.  An easy to use gui for simulating big data using tennessee eastman process.  _Quality and Reliability Engineering International_ , 38:264–282, 2 2022.  ISSN 0748-8017.  doi: 10.1002/qre.2975.
  * Burbidge et al. (1997) Robert Burbidge, Jem J Rowland, and Ross D King.  Active learning for regression based on query by committee.  In _Intelligent Data Engineering and Automated Learning - IDEAL 2007_ , 1997.
  * Cacciarelli and Kulahci (2022) Davide Cacciarelli and Murat Kulahci.  A novel fault detection and diagnosis approach based on orthogonal autoencoders.  _Computers & Chemical Engineering_, 2022.  doi: 10.1016/j.compchemeng.2022.107853.  URL https://www.sciencedirect.com/science/article/pii/S0098135422001910.
  * Cacciarelli and Kulahci (2023) Davide Cacciarelli and Murat Kulahci.  A survey on online active learning.  _arXiv preprint arXiv:2302.08893_ , 2023.  doi: 10.48550/arXiv.2302.08893.  URL https://arxiv.org/abs/2302.08893.
  * Cacciarelli et al. (2022) Davide Cacciarelli, Murat Kulahci, and John Sølve Tyssedal.  Stream-based active learning with linear models.  _Knowledge-Based Systems_ , 2022.  doi: 10.1016/j.knosys.2022.109664.  URL https://www.sciencedirect.com/science/article/pii/S0950705122008425.
  * Cacciarelli et al. (2023) Davide Cacciarelli, Murat Kulahci, and John Sølve Tyssedal.  Robust online active learning.  _arXiv preprint arXiv:2302.00422_ , 2023.  doi: 10.48550/arXiv.2302.00422.  URL https://arxiv.org/abs/2302.00422.
  * Cai et al. (2013) Wenbin Cai, Ya Zhang, and Jun Zhou.  Maximizing expected model change for active learning in regression.  _Proceedings - IEEE International Conference on Data Mining, ICDM_ , pages 51–60, 2013.  ISSN 15504786.  doi: 10.1109/ICDM.2013.104.
  * Capaci et al. (2019) Francesca Capaci, Erik Vanhatalo, Murat Kulahci, and Bjarne Bergquist.  The revised tennessee eastman process simulator as testbed for spc and doe methods.  _Quality Engineering_ , 31, 4 2019.  ISSN 0898-2112.  doi: 10.1080/08982112.2018.1461905.
  * de Almeida Moreira et al. (2021) Bruno Rafael de Almeida Moreira, Victor Hugo Cruz, Matheus Luís Cunha Oliveira, and Ronaldo da Silva Viana.  Full-scale production of high-quality wood pellets assisted by multivariate statistical process control.  _Biomass and Bioenergy_ , 151, 8 2021.  ISSN 18732909.  doi: 10.1016/j.biombioe.2021.106159.
  * Fortuna et al. (2007) Luigi Fortuna, Salvatore Graziani, Alessandro Rizzo, and Maria G Xibilia.  _Soft sensors for monitoring and control of industrial processes_ , volume 22.  Springer, 2007.
  * Frumosu and Kulahci (2018) Flavia D. Frumosu and Murat Kulahci.  Big data analytics using semi-supervised learning methods.  _Quality and Reliability Engineering International_ , 34:1413–1423, 11 2018.  ISSN 10991638.  doi: 10.1002/qre.2338.
  * Ge et al. (2014) Zhiqiang Ge, Biao Huang, and Zhihuan Song.  Nonlinear semisupervised principal component regression for soft sensor modeling and its mixture form.  _Journal of Chemometrics_ , 28:793–804, 11 2014.  ISSN 1099128X.  doi: 10.1002/cem.2638.
  * Grbić et al. (2013) Ratko Grbić, Dražen Slišković, and Petr Kadlec.  Adaptive soft sensor for online prediction and process monitoring based on a mixture of gaussian process models.  _Computers & Chemical Engineering_, 58:84–97, 11 2013\.  ISSN 00981354.  doi: 10.1016/j.compchemeng.2013.06.014.
  * Hotelling (1947) H Hotelling.  Multivariate quality control.  _Techniques of Statistical Analysis_ , 1947.  URL http://ci.nii.ac.jp/naid/10021322508/en/.
  * Jia et al. (2020) Xuqing Jia, Wende Tian, Chuankun Li, Xia Yang, Zhongjun Luo, and Hui Wang.  A dynamic active safe semi-supervised learning framework for fault identification in labeled expensive chemical processes.  _Processes_ , 8, 1 2020.  ISSN 22279717.  doi: 10.3390/pr8010105.
  * Jiang et al. (2017) Li Jiang, Zhiqiang Ge, and Zhihuan Song.  Semi-supervised fault classification based on dynamic sparse stacked auto-encoders model.  _Chemometrics and Intelligent Laboratory Systems_ , 168:72–83, 9 2017.  ISSN 18733239.  doi: 10.1016/j.chemolab.2017.06.010.
  * Kadlec et al. (2009) Petr Kadlec, Bogdan Gabrys, and Sibylle Strandt.  Data-driven soft sensors in the process industry.  _Computers & Chemical Engineering_, 33(4):795–814, 2009.  ISSN 0098-1354.  doi: https://doi.org/10.1016/j.compchemeng.2008.12.012.  URL https://www.sciencedirect.com/science/article/pii/S0098135409000076.
  * Reinartz et al. (2021) Christopher Reinartz, Murat Kulahci, and Ole Ravn.  An extended tennessee eastman simulation dataset for fault-detection and decision support systems.  _Computers & Chemical Engineering_, 149:107281, 6 2021\.  ISSN 00981354.  doi: 10.1016/j.compchemeng.2021.107281.
  * Ricker (1996) N. Lawrence Ricker.  Decentralized control of the tennessee eastman challenge process.  _Journal of Process Control_ , 6, 8 1996.  ISSN 09591524.  doi: 10.1016/0959-1524(96)00031-5.
  * Scott (1992) David W. Scott.  _Multivariate Density Estimation: Theory, Practice, and Visualization_.  John Wiley & Sons, Inc., 8 1992.  ISBN 9780471547709.  doi: 10.1002/9780470316849.
  * Settles (2009) Burr Settles.  Computer sciences department active learning literature survey, 2009.
  * Wang et al. (2019) Wei Wang, Dan Yang, Feiyu Chen, Yunsheng Pang, Sheng Huang, and Yongxin Ge.  Clustering with orthogonal autoencoder.  _IEEE Access_ , 7:62421–62432, 2019.  ISSN 21693536.  doi: 10.1109/ACCESS.2019.2916030.
  * Yuan et al. (2018) Xiaofeng Yuan, Biao Huang, Yalin Wang, Chunhua Yang, and Weihua Gui.  Deep learning-based feature representation and its application for soft sensor modeling with variable-wise weighted sae.  _IEEE Transactions on Industrial Informatics_ , 14:3235–3243, 7 2018.  ISSN 15513203.  doi: 10.1109/TII.2018.2809730.
  * Zhu et al. (2015) Jinlin Zhu, Zhiqiang Ge, and Zhihuan Song.  Robust semi-supervised mixture probabilistic principal component regression model development and application to soft sensors.  _Journal of Process Control_ , 32:25–37, 5 2015.  ISSN 09591524.  doi: 10.1016/j.jprocont.2015.04.015.


# Diamond quantum magnetometer with dc sensitivity of ¡ 10 pT Hz-1/2 toward
measurement of biomagnetic field

N. Sekiguchi<EMAIL_ADDRESS>Department of Electrical and
Electronic Engineering, Tokyo Institute of Technology, Meguro, Tokyo 152-8550,
Japan  M. Fushimi  Department of Bioengineering, The University of Tokyo,
Bunkyo, Tokyo 113-8656, Japan  A. Yoshimura  Department of Electrical and
Electronic Engineering, Tokyo Institute of Technology, Meguro, Tokyo 152-8550,
Japan  C. Shinei  Research Center for Electronic and Optical Materials,
National Institute for Materials Science, Tsukuba, Ibaraki 305-0044, Japan  M.
Miyakawa  T. Taniguchi  Research Center for Materials Nanoarchitectonics,
National Institute for Materials Science, Tsukuba, Ibaraki 305-0044, Japan  T.
Teraji  Research Center for Electronic and Optical Materials, National
Institute for Materials Science, Tsukuba, Ibaraki 305-0044, Japan  H. Abe  S.
Onoda  T. Ohshima  Takasaki Institute for Advanced Quantum Science, National
Institutes for Quantum Science and Technology, Takasaki, Gunma 370-1292, Japan
M. Hatano  Department of Electrical and Electronic Engineering, Tokyo
Institute of Technology, Meguro, Tokyo 152-8550, Japan  M. Sekino  Department
of Bioengineering, The University of Tokyo, Bunkyo, Tokyo 113-8656, Japan  T.
Iwasaki  Department of Electrical and Electronic Engineering, Tokyo Institute
of Technology, Meguro, Tokyo 152-8550, Japan

###### Abstract

We present a sensitive diamond quantum sensor with a magnetic field
sensitivity of $9.4\pm 0.1~{}\mathrm{pT/\sqrt{Hz}}$ in a near-dc frequency
range of 5 to 100 Hz. This sensor is based on the continuous-wave optically
detected magnetic resonance of an ensemble of nitrogen–vacancy centers along
the [111] direction in a diamond (111) single crystal. The long
$T_{2}^{\ast}\sim 2~{}\mathrm{\mu s}$ in our diamond and the reduced intensity
noise in laser-induced fluorescence result in remarkable sensitivity among
diamond quantum sensors. Based on an Allan deviation analysis, we demonstrate
that a sub-picotesla field of 0.3 pT is detectable by interrogating the
magnetic field for a few thousand seconds. The sensor head is compatible with
various practical applications and allows a minimum measurement distance of
about 1 mm from the sensing region. The proposed sensor facilitates the
practical application of diamond quantum sensors.

##  I  Introduction

The biomedical applications of quantum sensors have been studied for over a
decade [1]. The realization of magnetoencephalography (MEG) under ambient
conditions is a major goal (conventional MEG requires a magnetically shielded
room). In addition to clinical diagnosis [2, 3], ambient-condition MEG can be
used for daily diagnosis, brain-machine interfaces [4, 5], and fundamental
research on brain function [6, 7, 8, 9]. A quantum magnetometer that uses
nitrogen–vacancy (NV) centers in diamond is a candidate for realizing ambient-
condition MEG given that it can be operated with high sensitivity at room
temperature in an ambient magnetic field [10, 11, 12, 13, 14, 15, 16, 17, 18].
A spatial resolution on the millimeter scale or below, far better than the
centimeter-scale resolution of conventional MEG [3], is expected for a diamond
quantum magnetometer [19].

Magnetometry based on continuous-wave optically detected magnetic resonance
(CW-ODMR) is the most widely used method for measuring magnetic fields using
NV centers [10, 11, 13, 15, 14, 12, 16]. In this method, a microwave (MW)
field continuously drives the magnetic resonance of the NV center spin and the
spin state is continuously read out as the intensity of the laser-induced
fluorescence from the NV center. Compared with other methods based on pulsed
MWs and/or light [17, 13, 18, 20], the CW-ODMR method has a simpler
experimental setup and is easier to apply to actual measurements. Millimeter-
scale magnetocardiography [19] has been realized using the CW-ODMR method.
However, to realize MEG, measurement of an encephalomagnetic field requires
exceptional sensitivity (on the order of $\mathrm{pT/\sqrt{Hz}}$). The
frequency of a clinically relevant encephalomagnetic field ranges from nearly
dc to $\sim 100~{}{\mathrm{H}z}$ [3, 2, 9]. Reported sensitivities are worse
than required in this frequency range. For example, sensitivities of around
$20$ to $30~{}\mathrm{pT/\sqrt{Hz}}$ [13, 15] have been demonstrated. In
addition, a sensitivity of $15~{}\mathrm{pT\sqrt{Hz}}$ in a higher frequency
range (80 Hz to 3.6 kHz) has been reported [12]. A short standoff distance
from field generating sources in the brain is also required given that the
decay of an encephalomagnetic field is inversely proportional to the square of
the distance [2]. Therefore, for biomedical applications, the sensitivity in
the near-dc frequency range of a diamond quantum magnetometer that can closely
approach the target object must be improved.

Here, we develop a CW-ODMR-based diamond magnetometer for practical
applications (e.g., MEG of a living animal). The sensor head of the
magnetometer was designed to approach the target object to a distance of about
1 mm with a sensing volume of $0.03~{}\mathrm{mm^{3}}$. By carefully tuning
the experimental conditions and using a high-quality diamond, we achieved a
record-breaking sensitivity of $9.4\pm 0.1~{}\mathrm{pT/\sqrt{Hz}}$ in a near-
dc frequency range of 5 to 100 Hz. Based on the Allan deviation, the minimum
detectable field was found to be 8.5 and 0.3 pT for measurement periods of 1
second and several thousand seconds, respectively.

##  II  Experimental setup

###  II.1  Sensor head

In this work, we synthesized a single-crystalline diamond using a high-
pressure–high-temperature (HPHT) method with a 12C isotopically enriched
carbon source. The reduced concentration of 13C was about 500 ppm. The amount
of titanium in the metal solvent in the HPHT synthesis was adjusted to control
the initial concentration of neutral substituted nitrogen
($\mathrm{N_{s}^{0}}$) in the diamond crystal [21]. The initial
$[\mathrm{N_{s}^{0}}]$ was estimated to be 5.6 ppm using electron spin
resonance. The origin of nitrogen in diamond crystals seems to be impurities
introduced from the source material, solvent, or pressure transmitting medium
during the growth process. Since this nitrogen is of natural origin, the
isotope abundance of $\mathrm{N_{s}^{0}}$ is the same as the natural abundance
(14N, 99.6%; 15N, 0.4%). After this HPHT synthesis, a piece of the crystal was
cut out parallel to the (111) crystal plane. The dimensions of this diamond
sample were approximately $1~{}\mathrm{mm}\times 0.7~{}\mathrm{mm}$ in area
and 0.4 mm in thickness. Negatively charged NV ($\mathrm{NV^{-}}$) centers
were then produced using electron beam irradiation followed by annealing at
$1000~{}\mathrm{C^{\circ}}$ for 2 hours in vacuum. The energy and total
fluence of the irradiation were 2.0 MeV and $5\times
10^{17}~{}\mathrm{cm^{-2}}$, respectively. The concentrations of the produced
$\mathrm{NV^{-}}$ and residual $\mathrm{N_{s}^{0}}$ were estimated to be 1.2
and 2.3 ppm, respectively, using electron spin resonance [22]. A full width at
half maximum of 0.19 MHz for the CW-ODMR peak was experimentally measured
independent of this work. This linewidth indicates a long dephasing time of
$T_{2}^{\ast}\sim 2~{}\mathrm{\mu s}$.

Figure 1:  Experimental setup (not to scale). (a) Sensor head design and MW
circuit diagram. $\mathrm{PD_{fl}}$: fluorescence photodiode; LPF: long-pass
filter. (b) Optical setup. $\lambda/2$: halfwave plate; PBS: polarizing beam
splitter; NPBS: non-polarizing beam splitter; L: lens; M: mirror;
$\mathrm{PD_{ref}}$: reference photodiode; BB: beam block.

The conceptual design of our sensor head is shown in Fig. 1(a). This sensor
head was designed to closely approach the head of a living animal and measure
the encephalomagnetic field along the $z$ axis by an ensemble of NV centers
oriented to the surface-normal [111] direction parallel to the $z$ axis. The
sensor head components described in this section (see below) were integrated
using plastic and aluminum holders. Hence, the sensor head can be freely moved
as a unit and easily positioned close to the target object.

The diamond containing NV centers was attached by a high-thermal-conductivity
glue to a polycrystalline diamond plate ($10\times 10\times
0.5~{}\mathrm{mm^{3}}$) in order to dissipate the heat due to laser
illumination. The other side of the polycrystalline diamond plate had a
current flow guide for MWs. The MW guide was made of thin copper film; the
distance between the lower side of the MW guide and the excited NV centers was
0.8 mm. A bias magnetic field of 0.9 mT along the $z$ axis was applied by a
ring samarium-cobalt magnet.

We used a hemispherical lens with a high refractive index of 2.0 to enhance
the collection efficiency of the laser-induced fluorescence from the NV center
ensemble [14]. The fluorescence collection efficiency from the diamond surface
facing the lens [top surface in Fig. 1(a)] was assumed to be as high as about
56% based on a previously reported numerical calculation [14] for a similar
setup. The fluorescence that was not emitted from this surface was considered
to be emitted mainly from the side faces due to the high refractive index
(2.4) of the diamond [23]. Some of the fluorescence from the side faces of the
diamond was collected by the lens since the lens diameter (4 mm) was larger
than the size of the diamond. Fluorescence was also collected by an
elliptically shaped reflective inner surface of an aluminum block. Stray green
light and part of the fluorescence from neutrally charged NV
($\mathrm{NV^{0}}$) centers were filtered out by a long-pass filter with a
cut-on wavelength of 633 nm. The transmitted fluorescence was detected by a
reverse-biased photodiode ($\mathrm{PD_{fl}}$).

###  II.2  CW-ODMR measurement setup

The NV ensemble was excited by a green laser at 532 nm from a side face of the
diamond, as shown in Fig. 1(a). Figure 1(b) shows the optical setup. A laser
beam with a diameter of about 3 mm was focused by a lens with a focal length
of 300 mm. The beam diameter at the diamond was estimated to be
$200~{}\mathrm{\mu m}$. The excitation volume was estimated to be
$0.03~{}\mathrm{mm^{3}}$. The laser light was linearly polarized along the $y$
axis, which is perpendicular to the chosen NV orientation. The fluorescence
photocurrent $I_{\mathrm{fl}}=6.6~{}\mathrm{mA}$ was observed at an incident
light power of 0.39 W, which corresponds to a detected fluorescence power of
about 13 mW.

The noise in the fluorescence due to the intensity fluctuation of the incident
laser was reduced using a balanced detection technique. The reference light,
which was picked up by a non-polarizing beam splitter, was detected by a
reverse-biased photodiode ($\mathrm{PD_{ref}}$). In this work, we connected
the anode of $\mathrm{PD_{fl}}$ to the cathode of $\mathrm{PD_{ref}}$ to
obtain the difference between their photocurrents, $I_{\mathrm{fl}}$ and
$I_{\mathrm{ref}}$, respectively. The difference photocurrent
$I_{\mathrm{diff}}$ was amplified by a lab-built transimpedance amplifier with
a gain of 10 kV/A. The power of the reference light was finely adjusted using
a halfwave plate and a polarizing beam splitter to achieve a high reduction
rate for the intensity noise. The polarization fluctuation of the laser was
converted into an intensity fluctuation by a polarization beam splitter just
after the laser. The beam diameter at $\mathrm{PD_{ref}}$ was expanded by a
lens to balance the nonlinear response of the photodiode with that of
$\mathrm{PD_{fl}}$, since the nonlinear response depends on spot size [24].

The magnetic resonance between the ground states $|0\rangle$ and $|-1\rangle$
was driven by applying an MW current to the MW guide. To enhance the amplitude
of a CW-ODMR peak, we simultaneously drove the three transitions associated
with the hyperfine spin state by three-tone MWs [12], which was generated by
mixing radio-frequency (RF) waves at 2.16 MHz with MWs and summing the mixed
waves with bypassed MWs. In this work, the enhancement factor for the peak
amplitude was about 2.5.

We adopted lock-in detection, achieved by sinusoidally modulating the MW
frequency, to avoid large residual noise at low frequencies. The amplified
difference photocurrent was fed into a lock-in amplifier and demodulated with
the modulation frequency as the reference. The 3-dB cutoff frequency of the
low-pass filter in the lock-in amplifier was 149.4 Hz, which corresponds to a
noise-equivalent-power bandwidth of 168.8 Hz. The demodulated output was
recorded on a computer via an analog-to-digital converter.

The sensor head and optical setup were inside a room that was shielded from
magnetic fields by three permalloy layers to reduce environmental field
fluctuations. The total shielding factor of this room was about $2\times
10^{-4}$ at 1 Hz and about $1\times 10^{-5}$ at 10 Hz. An additional permalloy
box was placed around the sensor head. The front face of the shield box
remained open to introduce the incident laser and the target object.

##  III  Results

###  III.1  CW-ODMR measurement

Figure 2:  Lock-in CW-ODMR spectrum (a) over hyperfine manifold and (b) in
near-resonant region of central peak. Measured demodulated photocurrent is
shown by filled circles. The solid line in (a) represents the fitted curve
obtained from the summation of five derivative Lorentzian functions. The
linear function shown by the dashed line in (b) was fitted to the near-
resonant data to obtain a zero-crossing slope.

Figure 2(a) shows a CW-ODMR spectrum of the ensemble of [111]-oriented NV
centers. The vertical axis in the figure is the demodulated signal $\tilde{I}$
in the photocurrent, which was calculated using the gains at the
transimpedance and lock-in amplifiers. The horizontal axis is the detuning
$\delta$ from the resonance frequency of the central peak at which the three
hyperfine spin states were simultaneously driven. The fluorescence
photocurrent $I_{\mathrm{fl}}$ at a far-detuned MW frequency was
$I_{\mathrm{fl}}=6.6~{}\mathrm{mA}$. In this measurement, the frequency and
depth of the modulation were 6.2 kHz and 160 kHz, respectively. We found that
a modulation frequency of 3 to 7 kHz yielded a low-noise output. The
modulation frequency was finely tuned within this range on each day of the
experiment because the frequencies of some noise peaks in the intensity noise
slightly shifted over time. The modulation frequency of higher than several
kilohertz caused a decrease in the amplitude of the CW-ODMR peaks, probably
because the dynamics of the population in the ground states was slower than
the modulation at several kilohertz for the relatively low power of $<1$ W and
large spot size of 200 $\mathrm{\mu m}$ in diameter. The MW and RF wave power
was tuned to yield a maximum zero-crossing slope at the central peak. The
black solid curve is fitted to the measured data using the summation of five
derivative Lorentzian functions. The corresponding full width at half maximum
of the derivative Lorentzian function was about 0.48 MHz . This linewidth is
greater than the inhomogeneous broadening of the CW-ODMR peak due to the
inhomogeneity in the bias magnetic field, which was estimated to be
approximately 0.2 MHz.

To determine the zero-crossing slope, we measured a CW-ODMR spectrum at the
near-resonance region of the central peak, as shown in Fig. 2(b). The
demodulated photocurrent linearly depends on the detuning in this region. The
zero-crossing slope $d\tilde{I}/d\delta$ was measured to be
$324~{}\mathrm{pA/Hz}$ by fitting the data with a linear function, as shown by
the black dashed line. This slope corresponds to the photocurrent response to
magnetic field variation as
$(d\tilde{I}/d\delta)/\gamma_{\mathrm{e}}=9.06~{}\mathrm{A/T}$, where
$\gamma_{\mathrm{e}}=28.0~{}\mathrm{GHz/T}$ is the gyromagnetic ratio for an
NV center.

###  III.2  Reduction in intensity noise

The reduction rate for the intensity noise was estimated at
$I_{\mathrm{fl}}=25~{}\mathrm{mA}$. We measured the standard deviations of
$\tilde{I}$ with and without the reference light to be $3.0~{}\mathrm{nA}$ and
$130~{}\mathrm{nA}$, respectively. Here, the MW source was switched off to
isolate the sensor from the noise associated with the environmental magnetic
field. The relative intensity noise in the incident light was roughly
estimated to be
$10\mathrm{log_{10}}\left(\frac{130~{}\mathrm{nA}^{2}}{168.8~{}\mathrm{Hz}\times
25~{}\mathrm{mA}^{2}}\right)=-130~{}\mathrm{dBc/Hz}$ at a modulation frequency
of 6.2 kHz. The photon shot noise with and without the reference light was
calculated to be $1.6~{}\mathrm{nA}$ and $1.2~{}\mathrm{nA}$, respectively.
The details of the photon shot noise calculation are described in Sec. III.3.
We obtained the following reduction rate for the fluorescence intensity noise:

$\sqrt{\frac{3.0~{}\mathrm{nA}^{2}-1.6~{}\mathrm{nA}^{2}}{130~{}\mathrm{nA}^{2}-1.2~{}\mathrm{nA}^{2}}}=1.9\times
10^{-2}.$

We found that this “red–green” balance detection exhibited a similar reduction
rate to that for the “green–green” balanced detection with the reference and
incident lights.

###  III.3  Photocurrent dependence of noise

Figure 3:  Fluorescence photocurrent dependence of (a) floor of noise spectral
density of $\tilde{I}$, (b) zero-crossing slope $d\tilde{I}/d\delta$, and (c)
estimated floor of equivalent magnetic field noise spectral density.

We analyzed the noise components (photon shot noise, fluorescence intensity
noise, and electrical noise) of the detectors and circuits by measuring their
dependence on the fluorescence photocurrent $I_{\mathrm{fl}}$. The demodulated
photocurrent $\tilde{I}$ was recorded for 5 s and Fourier-transformed to
provide a single-sided noise amplitude spectral density $n_{\tilde{I}}$. To
evaluate the noise $n_{\tilde{I},\mathrm{far}}$ without influence from
environmental magnetic field noise, the analysis was performed with an MW
carrier frequency of 2.4 GHz, which was far-detuned from the resonance. We
observed no excess noise due to the application of the far-detuned MWs. The
noise density $n_{\tilde{I},\mathrm{far}}$ was almost flat up to the cutoff
frequency of the lock-in amplifier. The average $\langle
n_{\tilde{I},\mathrm{far}}\rangle$ of the noise density within the 100-Hz
bandwidth was taken as a measure of the intrinsic noise of our diamond sensor
at a given $I_{\mathrm{fl}}$. The dependence on $I_{\mathrm{fl}}$ of $\langle
n_{\tilde{I},\mathrm{far}}\rangle$ is shown in Fig. 3(a). Here, we varied the
incident laser power using a halfwave plate and a polarizing beam splitter
just before the non-polarizing beam splitter. In the figure, the measured
$\langle n_{\tilde{I},\mathrm{far}}\rangle$ is represented by open circles.
The relative uncertainty in the data, shown as error bars, was independently
evaluated to be 5%. $\langle n_{\tilde{I},\mathrm{far}}\rangle$ at
$I_{\mathrm{fl}}=0$ represents the electrical noise density $\langle
n_{\tilde{I},\mathrm{elec}}\rangle$ and was measured to be
$20~{}\mathrm{pA/\sqrt{Hz}}$ by blocking the laser beam before the non-
polarizing beam splitter.

We fitted the noise model in Eq. (1) to the data.

$\langle n_{\tilde{I},\mathrm{far}}\rangle=\sqrt{\langle
n_{\tilde{I},\mathrm{elec}}\rangle^{2}+p_{1}I_{\mathrm{fl}}+p_{2}I_{\mathrm{fl}}^{2}}.$
(1)

The second and third terms represent the photon shot noise $\langle
n_{\tilde{I},\mathrm{psn}}\rangle$ and fluorescence intensity noise $\langle
n_{\tilde{I},\mathrm{int}}\rangle$, respectively. This noise model well
describes the data, as shown by the black solid curve in Fig. 3(a). The fitted
parameters were $p_{1}=(5.0\pm 0.6)\times 10^{-19}~{}\mathrm{A/Hz}$ and
$p_{2}=(5.0\pm 0.5)\times 10^{-17}~{}\mathrm{/Hz}$. The black dashed curve is
the sum of $\langle n_{\tilde{I},\mathrm{elec}}\rangle$ and the calculated
shot noise given by

$\sqrt{\langle n_{\tilde{I},\mathrm{elec}}\rangle^{2}+2\times
2q_{e}I_{\mathrm{fl}}},$ (2)

where $q_{e}=1.6\times 10^{-19}~{}\mathrm{C}$ is the elementary charge. The
factor of 2 for the shot noise term was introduced because the shot noise at
the two photodiodes was assumed to be independent. The measured shot noise
coefficient $p_{1}=(5.0\pm 0.6)\times 10^{-19}~{}\mathrm{A/Hz}$ is close to
the calculated value of $2\times 2q_{e}=6.4\times 10^{-19}~{}\mathrm{A/Hz}$.
The intensity noise $\langle
n_{\tilde{I},\mathrm{int}}\rangle=\sqrt{p_{2}I_{\mathrm{fl}}^{2}}$ is
equivalent to $\langle n_{\tilde{I},\mathrm{psn}}\rangle$ at the fluorescence
photocurrent $I_{\mathrm{fl,eqv}}=p_{1}/p_{2}=10\pm 1.6~{}\mathrm{mA}$. The
photon shot noise surpassed the laser intensity noise at a low fluorescence
photocurrent ($<I_{\mathrm{fl,eqv}}$).

The sensor noise $n_{B}$ in the magnetic field measurement depends on
demodulated photocurrent noise $n_{\tilde{I}}$ and zero-crossing slope
$d\tilde{I}/d\delta$ as $n_{B}=n_{\tilde{I}}/(\gamma_{e}d\tilde{I}/d\delta)$.
The fluorescence dependence of the slope was measured, as shown in Fig. 3(b).
The error bars are the estimated standard deviations of the slope; they are
much smaller than the marker size. Here, the modulation parameters and powers
of the MWs and RF waves were fixed over all measurements; they were tuned at
$I_{\mathrm{fl}}=7.2~{}\mathrm{mA}$ to maximize the slope. Note that the
optimal parameters and power depend on the incident laser power [25, 26].
Nevertheless, we confirmed that tuning these values resulted in an improvement
in the slope of about 3% at $I_{\mathrm{fl}}=29.3~{}\mathrm{mA}$. We thus
assumed that the relative uncertainty of the data was several percent in this
measurement.

We found that the slope saturated as the incident laser power increased. This
saturation could be explained by a charge state conversion of NV centers,
which led to a decrease in the contrast of a CW-ODMR peak, because
[$\mathrm{N_{s}^{0}}$] was only about twice as large as [$\mathrm{NV^{-}}$] in
our diamond [27, 28, 29]. A detailed investigation of this saturation is
beyond the scope of this work. The magnetic field noise density $\langle
n_{B,\mathrm{far}}\rangle$ expected from the measured $\langle
n_{\tilde{I},\mathrm{far}}\rangle$ and $d\tilde{I}/d\delta$ did not
monotonically decrease as $I_{\mathrm{fl}}$ increased, as shown in Fig. 3(c),
because of the saturation of the slope. The error bars indicate the
uncertainties computed from a relative uncertainty of 5% in the slope and the
covariance matrix used in the curve fitting to $\langle
n_{\tilde{I},\mathrm{far}}\rangle$ with Eq. (1). The photocurrent dependence
of $\langle n_{B,\mathrm{far}}\rangle$ suggests that good sensitivity to a
magnetic field can be achieved at $I_{\mathrm{fl}}$ from 5 to 20 mA.

###  III.4  Magnetic field noise spectral density and sensitivity

Figure 4:  Single-sided noise amplitude spectral density in magnetic field
measurement. The blue and orange curves are the sensor noise measured with
resonant and far-detuned MWs, respectively. Calculated photon-shot-noise-
limited sensitivity of $6.9~{}\mathrm{pT/\sqrt{Hz}}$ is indicated by the
dashed line. PSN: photon shot noise.

We measured single-sided noise amplitude spectral density $n_{B,\mathrm{res}}$
in a magnetic field measurement where the MWs were resonant with the central
CW-ODMR peak ($\delta=0$). The noise spectrum $n_{B,\mathrm{res}}$ was
computed using the discrete Fourier transform from a measured time trace of
$\tilde{I}$ for 5 s with sampling frequency $F_{s}=400~{}\mathrm{Hz}$. Figure
4 shows the measured $n_{B,\mathrm{res}}$, which was averaged over 10 time
measurements, at $I_{\mathrm{fl}}=6.4~{}\mathrm{mA}$ (blue solid curve). The
optimal power of the reference light was estimated from the CW-ODMR peak
contrast of 3% and a reference light power that had been optimized with the
far-detuned MWs. The zero-crossing slope $d\tilde{I}/d\delta$ was $332\pm
0.7~{}\mathrm{pA/Hz}$. Note that the displayed $n_{B,\mathrm{res}}$ was
digitally filtered by narrow-band notch filters for the harmonics of a 50-Hz
power line and a band-pass filter with 3-dB cutoff frequencies of 5 and 100
Hz, which corresponds to the bandwidth of the target object (e.g., brain of a
living animal). The noise-equivalent power bandwidth $f_{\mathrm{NEP}}$ of the
digital filtering was numerically calculated to be 91.9 Hz. In this numerical
calculation, white noise with a standard deviation of $\sigma$ was numerically
computed with sampling frequency $F_{s}$ and digitally filtered. The standard
deviation $\sigma^{\prime}$ of the filtered noise, given by
$\sigma^{\prime}=\sigma\sqrt{f_{\mathrm{NEP}}/(F_{s}/2)}$ [30], was
numerically calculated to yield $f_{\mathrm{NEP}}$.

The achieved noise density indicated a very low floor in the single-sided
spectrum in the near-dc range. The lowest noise density floor, about
$9~{}\mathrm{pT/\sqrt{Hz}}$, was measured near 40 Hz and from 70 to 90 Hz. The
sudden drop in $n_{B,\mathrm{res}}$ at 90 Hz was due to the digital band-pass
filter. A low noise density of $15~{}\mathrm{pT/\sqrt{Hz}}$ was obtained even
near 5 Hz, even though magnetic field noise generally deteriorates at lower
frequency [13, 15, 19, 12, 20, 16]. We attributed the noise peaks around 25 Hz
to the mechanical vibration of the sensor head since these peaks shifted to
lower frequencies when the sensor head was additionally supported. The noise
spectral density $n_{B,\mathrm{far}}$ measured with the far-detuned MWs at
$I_{\mathrm{fl}}=6.4~{}\mathrm{mA}$ is shown by the orange trace in Fig. 4. We
found that $n_{B,\mathrm{far}}$ reached the photon-shot-noise-limited
sensitivity of $6.9~{}\mathrm{pT/\sqrt{Hz}}$ (black dashed line).

The sensitivity of our sensor was evaluated from the average power of the
measured noise, which was digitally filtered. The sensitivity $\eta$ is
defined as $\eta=\delta B\sqrt{T}$, where $\delta B$ is the minimum detectable
magnetic field for measurement time $T$. In this definition, the noise
spectrum is assumed to be frequency-independent (white noise). The standard
deviation of $128\pm 2~{}\mathrm{pT}$ in the measured time trace is considered
to represent $\delta B$ for measurement time $T=F_{s}^{-1}=2.5~{}\mathrm{ms}$.
Since the bandwidth of the digital band-pass filter was narrower than the
measurement bandwidth $F_{s}/2$ and the lock-in amplifier’s bandwidth,
$f_{\mathrm{NEP}}$ for the digital filtering was substituted for the
measurement bandwidth; that is, the sensitivity was equivalent to $\eta=\delta
B/\sqrt{2f_{\mathrm{NEP}}}$ [12, 30]. We achieved a sensitivity of
$\eta=9.4\pm 0.1~{}\mathrm{pT/\sqrt{Hz}}$.

###  III.5  Allan deviation

Figure 5:  Allan deviation as function of averaging time. Open circles show
calculated overlapping Allan deviation from a continuous measurement for 200
minutes. Estimated uncertainties of the Allan deviations are indicated by the
error bars, which are much smaller than the marker circle size.

The Allan deviation of the noise for a measurement of about 200 minutes was
computed to evaluate the stability of our sensor. We continuously tuned the MW
carrier frequency to the resonance using the demodulated photocurrent
$\tilde{I}$ output, which was low-pass-filtered with a cutoff frequency of 10
Hz. The bandwidth of the feedback response was approximately 2 Hz. The
measured noise was recorded every minute on a computer. The notch and band-
pass digital filters used in the sensitivity analysis (see Sec. III.4) were
not used in this analysis. We then computed the overlapping Allan deviation,
as shown in Fig. 5. The open circles indicate Allan deviations for a given
averaging time. The error bars represent the standard deviations of the Allan
deviations; they are much smaller than the marker size.

We found that the 1-second interrogation yielded an Allan deviation of 8.5 pT,
which is consistent with the evaluated sensitivity
($\eta=9.4~{}\mathrm{pT/\sqrt{Hz}}$). The Allan deviation showed a bump around
the 10-second averaging time. This bump may arise from a periodic fluctuation
of several tens of seconds. We found that the bump could be suppressed to some
extent by removing a 0.025-Hz component using a digital notch filter. However,
the cause of this fluctuation was not identified. The dashed line shows the
minimum detectable magnetic field $\delta B$ predicted by $\delta
B=\eta/\sqrt{T}$ with the achieved sensitivity of
$\eta=9.4~{}\mathrm{pT/\sqrt{Hz}}$. It indicates that the Allan deviation
scaled to $T^{-1/2}$ and was close to $\delta B$ at averaging times of 100 to
1000 seconds. The Allan deviation reached 0.3 pT for an averaging time of a
few thousand seconds and then seemed to saturate. The zero-crossing slopes
before and after this Allan deviation measurement were found to be almost the
same. We thus conclude that our sensor remained stable and could measure a
magnetic field with a sensitivity of $\eta=9.4\pm 0.1~{}\mathrm{pT/\sqrt{Hz}}$
for at least 200 minutes.

##  IV  Discussion

The demonstrated sensitivity of $9.4\pm 0.1~{}\mathrm{pT/\sqrt{Hz}}$ is the
best reported value for diamond quantum sensors based on the CW-ODMR of an
ensemble of NV centers [12, 30, 13, 15] in the frequency range of 5 to 100 Hz.
The previous best sensitivities were around $20$ to
$30~{}\mathrm{pT/\sqrt{Hz}}$ in the low frequency range [13, 15] and
$15~{}\mathrm{pT/\sqrt{Hz}}$ in the relatively high frequency range of 80
Hz–3.6 kHz [12]. Moreover, in our study, the noise floor of $n_{B}$ stayed
below $20~{}\mathrm{pT/\sqrt{Hz}}$ even at 5 Hz, as shown in Fig. 4. Given
that the noise environment is generally cleaner at higher frequency, the very
low noise floor of about $9~{}\mathrm{pT/\sqrt{Hz}}$ will continue into the
kilohertz range if we use a higher cutoff frequency for the lock-in
amplifier’s low-pass filter. The Allan deviation analysis showed that our
diamond magnetometer can interrogate a magnetic field for a long time with
remarkable sensitivity. Therefore, our sensor is capable of detecting a
repetitive biomagnetic field, for example, a stimulus-evoked field, with a
strength on the order of 1 pT by accumulating the signals.

CW-ODMR-based magnetometry has advantages over pulsed-MW-based magnetometry
for practical applications; it has a simpler experimental setup and looser
requirements for the inhomogeneities of the bias magnetic field and MWs.
Additionally, the use of a single orientation of NV center axes in our
magnetometer leads to a lower requirement for the bias field alignment
compared with that for multiple orientations [12, 15, 20]. The sensor head
design, which can approach the target object to a distance of 1 mm, relies on
a simple setup and reduced requirements for the bias field. The simplified
geometry between the single orientation of NV centers and the magnetic field
to be measured also facilitates various practical applications. We note that a
better sensitivity of around $2~{}\mathrm{pT/\sqrt{Hz}}$ in the low-frequency
range (from 10 Hz), achieved using the Ramsey method, has been recently
reported [20]; however, our sensor is more suitable for practical applications
such as biomagnetic field measurement because of its simplified setup and
short measurement distance.

We attributed a major part of the sensitivity improvement in this work to the
long dephasing time of $T_{2}^{\ast}\sim 2~{}\mathrm{\mu s}$, achieved by
decreasing the concentration of 13C to about 500 ppm and using a relatively
low initial nitrogen concentration of 5.6 ppm. The narrow linewidth of a CW-
ODMR peak due to the long $T_{2}^{\ast}$ resulted in a high response signal to
a magnetic field of $\gamma_{e}(d\tilde{I}/d\delta)=9.3~{}\mathrm{A/T}$, even
with the use of only a single crystallographic orientation of NV centers. The
photon-shot-noise-limited sensitivity is comparable to previously reported
values [12, 15]. In addition, the approximately five-fold improvement in the
intensity noise reduction in our balanced detection over the balanced
detection reported in a previous study [15] contributed to the good
sensitivity.

A reduction in the relative intensity noise of a laser can enhance
sensitivity. The relative intensity noise of our laser (Coherent Verdi G5) was
measured to be about $-130~{}\mathrm{dBc/Hz}$ at a modulation frequency of 6.2
kHz; the typical estimated relative intensity noise for state-of-the-art
solid-state lasers at the same frequency is $-140~{}\mathrm{dBc/Hz}$ [31, 32].
Therefore, a 10-dB improvement in $n_{\tilde{I},\mathrm{int}}^{2}$ is
feasible. This would result in a photon-shot-noise-limited sensitivity at up
to $I_{\mathrm{fl}}\sim 100~{}\mathrm{mA}$.

It is expected that the zero-crossing slope can be increased by extending the
dephasing time $T_{2}^{\ast}$ for the diamond. For example, a very long
dephasing time of $8.5~{}\mathrm{\mu s}$ with
$\mathrm{[NV^{-}]}=0.4~{}\mathrm{ppm}$ has been reported [13]. This long
dephasing time will offer a four-fold improvement if the same fluorescence
intensity is available since the shot-noise-limited sensitivity is
proportional to the linewidth of a CW-ODMR peak [25, 12]. Although the lower
$[\mathrm{NV^{-}}]$ emits weaker fluorescence, a photocurrent of up to 10 mA
can be obtained by increasing the incident laser power. In addition, the
fluorescence collection efficiency can be boosted to approximately unity by
using a total internal reflection lens and a light pipe [33, 20].

##  V  Conclusions

We demonstrated a sensitive diamond magnetometer with a magnetic field
sensitivity of $9.4\pm 0.1~{}\mathrm{pT/\sqrt{Hz}}$ in a near-dc frequency
range of 5 to 100 Hz. The magnetometer can closely approach the target object
and the measurement distance from the sensing volume was about 1 mm. The Allan
deviation indicated that our magnetometer can measure magnetic fields of 8.5
and 0.3 pT with a unity signal-to-noise ratio by interrogating for 1 second
and several thousands of seconds, respectively. Our high-sensitivity diamond
magnetometer was designed to be compatible with practical applications,
including the measurement of the encephalomagnetic field of a living animal.
The sensitivity improvement achieved in this work is an important step toward
realizing magnetoencephalography under ambient conditions with millimeter-
scale spatial resolution.

###### Acknowledgements.

This work was supported by the MEXT Quantum Leap Flagship Program (MEXT
Q-LEAP) Grant No. JPMXS0118067395 and JPMXS0118068379.

## References

  * Aslam  _et al._ [2023] N. Aslam, H. Zhou, E. K. Urbach, M. J. Turner, R. L. Walsworth, M. D. Lukin, and H. Park, Nature Reviews Physics 5, 157 (2023).
  * Hari  _et al._ [2018] R. Hari, S. Baillet, G. Barnes, R. Burgess, N. Forss, J. Gross, M. Hämäläinen, O. Jensen, R. Kakigi, F. Mauguière, N. Nakasato, A. Puce, G.-L. Romani, A. Schnitzler, and S. Taulu, Clinical Neurophysiology 129, 1720 (2018).
  * Uhlhaas  _et al._ [2018] P. J. Uhlhaas, T. Grent-’t-Jong, and J. Gross, JAMA Psychiatry 75, 969 (2018).
  * Fukuma  _et al._ [2018] R. Fukuma, T. Yanagisawa, H. Yokoi, M. Hirata, T. Yoshimine, Y. Saitoh, Y. Kamitani, and H. Kishima, Frontiers in Neuroscience 12, 478 (2018).
  * Rathee  _et al._ [2019] D. Rathee, A. Chowdhury, Y. K. Meena, A. Dutta, S. McDonough, and G. Prasad, IEEE Transactions on Neural Systems and Rehabilitation Engineering 27, 1020 (2019).
  * Boto  _et al._ [2018] E. Boto, N. Holmes, J. Leggett, G. Roberts, V. Shah, S. S. Meyer, L. D. Muñoz, K. J. Mullinger, T. M. Tierney, S. Bestmann, G. R. Barnes, R. Bowtell, and M. J. Brookes, Nature 555, 657 (2018).
  * Hill  _et al._ [2019] R. M. Hill, E. Boto, N. Holmes, C. Hartley, Z. A. Seedat, J. Leggett, G. Roberts, V. Shah, T. M. Tierney, M. W. Woolrich, C. J. Stagg, G. R. Barnes, R. Bowtell, R. Slater, and M. J. Brookes, Nature Communications 10, 4785 (2019).
  * Gross [2019] J. Gross, Neuron 104, 189 (2019).
  * Baillet [2017] S. Baillet, Nature Neuroscience 20, 327 (2017).
  * Taylor  _et al._ [2008] J. M. Taylor, P. Cappellaro, L. Childress, L. Jiang, D. Budker, P. R. Hemmer, A. Yacoby, R. Walsworth, and M. D. Lukin, Nature Physics 4, 810 (2008).
  * Acosta  _et al._ [2009] 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, Physical Review B 80, 115202 (2009).
  * Barry  _et al._ [2016] J. F. Barry, M. J. Turner, J. M. Schloss, D. R. Glenn, Y. Song, M. D. Lukin, H. Park, and R. L. Walsworth, Proceedings of the National Academy of Sciences 113, 14133 (2016).
  * Zhang  _et al._ [2021] C. Zhang, F. Shagieva, M. Widmann, M. Kübler, V. Vorobyov, P. Kapitanova, E. Nenasheva, R. Corkill, O. Rhrle, K. Nakamura, H. Sumiya, S. Onoda, J. Isoya, and J. Wrachtrup, Physical Review Applied 15, 064075 (2021).
  * Shim  _et al._ [2022] J. H. Shim, S.-J. Lee, S. Ghimire, J. I. Hwang, K.-G. Lee, K. Kim, M. J. Turner, C. A. Hart, R. L. Walsworth, and S. Oh, Physical Review Applied 17, 014009 (2022).
  * Graham  _et al._ [2023] S. Graham, A. Rahman, L. Munn, R. Patel, A. Newman, C. Stephen, G. Colston, A. Nikitin, A. Edmonds, D. Twitchen, M. Markham, and G. Morley, Physical Review Applied 19, 044042 (2023).
  * Fescenko  _et al._ [2020] I. Fescenko, A. Jarmola, I. Savukov, P. Kehayias, J. Smits, J. Damron, N. Ristoff, N. Mosavian, and V. M. Acosta, Physical Review Research 2, 023394 (2020).
  * Balasubramanian  _et al._ [2009] G. Balasubramanian, P. Neumann, D. Twitchen, M. Markham, R. Kolesov, N. Mizuochi, J. Isoya, J. Achard, J. Beck, J. Tissler, V. Jacques, P. R. Hemmer, F. Jelezko, and J. Wrachtrup, Nature Materials 8, 383 (2009).
  * Wolf  _et al._ [2015] T. Wolf, P. Neumann, K. Nakamura, H. Sumiya, T. Ohshima, J. Isoya, and J. Wrachtrup, Physical Review X 5, 041001 (2015).
  * Arai  _et al._ [2022] K. Arai, A. Kuwahata, D. Nishitani, I. Fujisaki, R. Matsuki, Y. Nishio, Z. Xin, X. Cao, Y. Hatano, S. Onoda, C. Shinei, M. Miyakawa, T. Taniguchi, M. Yamazaki, T. Teraji, T. Ohshima, M. Hatano, M. Sekino, and T. Iwasaki, Communications Physics 5, 200 (2022).
  * Barry  _et al._ [2023] J. F. Barry, M. H. Steinecker, S. T. Alsid, L. M. Pham, M. F. O’Keefe, and D. A. Braje, Sensitive AC and DC Magnetometry with Nitrogen-Vacancy Center Ensembles in Diamond (2023), arXiv:2305.06269 .
  * Miyakawa  _et al._ [2022] M. Miyakawa, C. Shinei, and T. Taniguchi, Japanese Journal of Applied Physics 61, 045507 (2022).
  * Shinei  _et al._ [2022] C. Shinei, Y. Masuyama, M. Miyakawa, H. Abe, S. Ishii, S. Saiki, S. Onoda, T. Taniguchi, T. Ohshima, and T. Teraji, Journal of Applied Physics 132, 214402 (2022).
  * Le Sage  _et al._ [2012] D. Le Sage, L. M. Pham, N. Bar-Gill, C. Belthangady, M. D. Lukin, A. Yacoby, and R. L. Walsworth, Physical Review B 85, 121202 (2012).
  * Scholze  _et al._ [2004] F. Scholze, R. M. Klein, and R. Mueller, Emerging lithographic technologies VIII 5374, 926 (2004).
  * Dréau  _et al._ [2011] A. Dréau, M. Lesik, L. Rondin, P. Spinicelli, O. Arcizet, J.-F. Roch, and V. Jacques, Physical Review B 84, 195204 (2011).
  * Jensen  _et al._ [2013] K. Jensen, V. M. Acosta, A. Jarmola, and D. Budker, Physical Review B 87, 014115 (2013).
  * Manson and Harrison [2005] N. Manson and J. Harrison, Diamond and Related Materials 14, 1705 (2005).
  * Collins [2002] A. T. Collins, Journal of Physics: Condensed Matter 14, 3743 (2002).
  * Lawson  _et al._ [1998] S. C. Lawson, D. Fisher, D. C. Hunt, and M. E. Newton, Journal of Physics: Condensed Matter 10, 6171 (1998).
  * Schloss  _et al._ [2018] J. M. Schloss, J. F. Barry, M. J. Turner, and R. L. Walsworth, Physical Review Applied 10, 034044 (2018).
  * Vernaleken  _et al._ [2012] A. Vernaleken, B. Schmidt, M. Wolferstetter, T. W. Hansch, R. Holzwarth, and P. Hommelhoff, Optics Express 20, 18387 (2012).
  * Vernaleken  _et al._ [2014] A. Vernaleken, B. Schmidt, T. W. Hänsch, R. Holzwarth, and P. Hommelhoff, Applied Physics B 117, 33 (2014).
  * Alsid  _et al._ [2023] S. T. Alsid, J. M. Schloss, M. H. Steinecker, J. F. Barry, A. C. Maccabe, G. Wang, P. Cappellaro, and D. A. Braje, Physical Review Applied 19, 054095 (2023).