# Hybrid Visual SLAM for Underwater Vehicle Manipulator Systems

Gideon Billings1, Richard Camilli2, Matthew Johnson-Roberson3 This work was
supported by NASA grant NNX16AL08G and by the National Science Foundation
under grants IIS-1830660 and IIS-1830500.1Gideon Billings is with Department
of Naval Architecture and Marine Engineering, University of Michigan, 2600
Draper Dr. Ann Arbor, MI 48109, USA<EMAIL_ADDRESS>Camilli is with
Applied Ocean Physics and Engineering, Woods Hole Oceanographic Institution,
Deep Submergence Laboratory Woods Hole, MA 02543, USA3Matthew Johnson-Roberson
is with The Robotics Institute, Carnegie Mellon University, 5000 Forbes Avenue
Pittsburgh, PA 15213-3890,
USA4https://github.com/gidobot/UWslam_datasetDigital Object Identifier (DOI):
10.1109/LRA.2022.3176448.

###### Abstract

This paper presents a novel visual feature based scene mapping method for
underwater vehicle manipulator systems (UVMSs), with specific emphasis on
robust mapping in natural seafloor environments. Our method uses GPU
accelerated SIFT features in a graph optimization framework to build a feature
map. The map scale is constrained by features from a vehicle mounted stereo
camera, and we exploit the dynamic positioning capability of the manipulator
system by fusing features from a wrist mounted fisheye camera into the map to
extend it beyond the limited viewpoint of the vehicle mounted cameras. Our
hybrid SLAM method is evaluated on challenging image sequences collected with
a UVMS in natural deep seafloor environments of the Costa Rican continental
shelf margin, and we also evaluate the stereo only mode on a shallow reef
survey dataset. Results on these datasets demonstrate the high accuracy of our
system and suitability for operating in diverse and natural seafloor
environments. We also contribute these datasets3 for public use.

###### Index Terms:

SLAM, Sensor Fusion, Computer Vision for Automation, Computer Vision for Other
Robotic Applications, Data Sets for Robotic Vision

##  I Introduction

Exploration vehicles for remote environments, such as rovers, planetary
landers, or underwater Remotely Operated Vehicles (ROVs) are often equipped
with manipulator systems for collecting samples, placing sensors, or otherwise
interacting with the environment. These systems largely rely on direct tele-
operation or manually scripted commands to execute manipulation tasks, due to
the risks associated with acting in unstructured and often complex remote
environments. Despite these risks, there are some remote environments, such as
Europa, the icy moon of Jupiter, so distant that any kind of tele-operation or
pre-scripted manipulator control is highly impractical. Considering
environments closer to home, the deep ocean is a domain of intensive
scientific research and exploration, but operation of deep submergence ROVs
with attendant support ships and pilots is costly and availability is limited.
Also, human error during tele-operated control can lead to damage to the
vehicle or loss of valuable mission time. These considerations motivate the
automation of manipulator systems for exploration vehicles, to enable complex
workspace interactions in communication limited or denied environments and to
support tele-operated procedures. This advancement can reduce the overhead
requirements of supporting tele-operation infrastructure (including associated
human labor and hotel functions), thereby increasing the availability of
robotic systems for scientific research. Critical to achieving safe and robust
autonomy of such vehicle-manipulator systems is scene perception and
reconstruction. In this work, we address the problem of feature based 3D scene
mapping for underwater vehicle-manipulator systems (UVMSs). A key innovation
of our mapping system is the fusion of feature points from both a vehicle
mounted stereo camera and a dynamically positioned manipulator mounted fisheye
camera into the same mapping framework. In situations where a UVMS’s movement
is limited or risky, our method addresses the challenge of limited viewpoints
from the vehicle mounted cameras and incomplete scene reconstruction due to
shadowing from scene structure by enabling the wrist mounted camera to
dynamically extend the map beyond the vehicle fixed camera views and fill in
shadowed areas of the scene.

This work makes the following contributions:

  1. 1.

To our knowledge, the first SLAM system, designed for manipulator systems,
that fuses a manipulator mounted fisheye camera into the same map with a
vehicle mounted stereo camera.

  2. 2.

An adaptation of the ORB-SLAM2 framework to GPU accelerated SIFT features,
with improved odometer based tracking and real-time performance.

  3. 3.

An evaluation of our method on both shallow reef and natural deep seafloor
environments, where our method achieves good performance and standard ORB-
SLAM2 fails. The evaluation datasets are also published with this paper4.

##  II Related Work

3D scene mapping is a very mature problem in computer vision and robotics, and
a rich body of literature has been generated from decades of study on the
topic. Here we present a review of the works which we consider most relevant
to our developed method and from which we took inspiration in our approach.

###  II-A Feature Based Visual SLAM

Since its inception, ORB-SLAM2 [1], remains one of the most widely adopted and
complete feature based SLAM systems, demonstrating that a bundle adjustment
approach can attain more accurate camera localization than direct methods or
ICP, with the advantage of being less computationally expensive. Given the
established robustness of ORB-SLAM2 across a variety of applications and
camera systems, the efficient computational performance based on a parallel
thread architecture, and the demonstrated accuracy of keyframe based bundle
adjustment for pose estimation, we chose to develop our method based on this
framework.

CoSLAM [2] proposed an innovative solution for fusing multiple synchronized
but independently moving monocular cameras into a single framework that can
also differentiate between dynamic and static feature points. We took
inspiration from this approach in our method design, with the key differences
being our use of stereo features to constrain the map scale, our fusion of
independent hybrid camera frames into the same map (i.e. the manipulator
mounted fisheye camera and a vehicle mounted perspective stereo camera), and
the specific adaptations of our method to underwater environments.

###  II-B Underwater SLAM

Significant progress has been made in underwater vision applied to large scale
survey reconstructions [3], terrain aided navigation [4, 5, 6], and ship hull
inspection [7]. However, dense scene reconstruction methods generally process
the image data offline, and methods designed for navigation generally provide
very sparse feature maps if any. In contrast, our method emphasises real-time
scene mapping, suitable for natural seafloor environments, that is robust to
underwater visual effects and provides an optimized feature map and camera
pose graph that can underlie dense reconstruction methods.

[4] proposed a stereo based SLAM method specifically designed for operating in
underwater feature-poor environments. The map is constructed as a pose graph
connecting to feature clusters. For inter-frame pose estimation of non
keyframes, they used the VISO2 stereo odometer [8], which they found to
perform better than the tracking stage in ORB-SLAM. For detecting loop
closures, this method generated a HALOC [9] signature for each feature
cluster, which can be efficiently matched across very large image sets and
does not require a prior training step like a bag of words representation.
This work informed our choice of using a modified version of VISO2 for the
initial inter-frame pose estimations. While their method was tailored
specifically to the problem of localization through the optimization of
keypoint cluster locations, our method, based on ORB-SLAM2, optimizes the
location of the individual map points, which is desirable for scene
reconstruction. [6] studied off the shelf monocular ORB-SLAM applied in
different shallow oceanic underwater environments. Their results showed that
ORB-SLAM performed well in structured or feature rich environments with
adequate lighting and low flickering. However, ORB-SLAM performed poorly in
areas with highly dynamic lighting, large numbers of moving objects, or low
textured regions such as sand beds.

###  II-C Kinematics in SLAM

Prior research has utilized eye-in-hand based SLAM, where a camera is mounted
near the end effector of a manipulator. ARM-SLAM [10] used a Kinect depth
sensor mounted on a manipulator with a fixed base to capture point clouds of
the scene and fuse them into a reconstruction using a method based on Kinect
Fusion. SKCLAM [11] used feature based pose tracking with an RGB-D camera on
the endeffector to calibrate the full kinematic parameters of an industrial
manipulator with a fixed base. Point clouds from the RGB-D camera were
integrated to construct a 3D map. [12] used ORB-SLAM3 and a stereo camera on a
mobile manipulator to map an orchard. Unlike these prior works, our method
fuses features from both an independent manipulator mounted fisheye camera and
a vehicle mounted stereo in a common feature graph. We use a monocular camera
on the wrist rather than relying on a depth sensor, which would be very bulky
to fit in a pressure rated housing for mounting on the manipulator. [13]
proposed a method for calibrating a dynamic camera cluster, where one camera
is articulated with respect to the other cameras in the system. They
demonstrated multi-camera SLAM with one camera mounted on a pan-tilt unit,
assuming accurate calibration of the pan-tilt unit’s extrinsics. In contrast,
our method is demonstrated with the manipulator camera having 5-DoF actuation
using very high baseline to the vehicle camera, and without relying on
accurate extrinsic measurement of the articulated camera. To our knowledge,
this is the first successful demonstration of an eye-in-hand SLAM method on
mobile underwater manipulator platforms in natural deep ocean environments.

##  III Method

In this section, we highlight the innovations made to adapt the ORB-SLAM2
system to SIFT features, the underwater environment, and the hybrid imagery.
For details on the system architecture that remain unchanged from ORB-SLAM2,
we defer the reader to [1].

Figure 1 shows a high level block diagram of the hybrid SLAM system, where our
method retains the same four threaded architecture as the original ORB-SLAM2.
Crucial modifications were made in the tracking thread, which follows the top
horizontal flow of the diagram, with separate functional flow branches for
stereo and monocular fisheye frames. Both stereo and fisheye frames share a
common keyframe representation which is processed through the local mapping,
loop closing, and full bundle adjustment threads. The core of our system is
the feature based stereo mapping framework, which can be operated stand-alone
or in a hybrid mode, where frames from an independently moving fisheye camera
are fused into the same map.

The constructed map is represented as a covisibilty graph of optimized
keyframe and keypoint poses, with factors between keyframes formed through
common keypoint observations. Like ORB-SLAM2, the covisibility graph is used
to retrieve a local neighborhood of keypoints for the tracking and local
mapping stages and forms the graph structure for the bundle adjustment
optimizations. A minimum spanning tree is also maintained, which connects
every keyframe to the neighbor with the maximum number of shared keypoint
observations. The spanning tree is used to propagate keyframe pose
optimizations from full bundle adjustment to new keyframes that were not
included during the optimization. A DBoW2 module [14], that we adapted to SIFT
features, is used for place recognition during relocalization and loop
closing.

Figure 1: System block diagram

###  III-A Hybrid Camera System

The hybrid camera system is specifically tailored to mobile manipulator
systems, with a stereo camera mounted on the vehicle and an independent but
synchronized fisheye camera mounted on the manipulator wrist. In our
evaluations, the stereo pair uses a pinhole camera model, calibrated with the
ROS [15] camera calibration package, and the fisheye camera uses the Kannala-
Brandt [16] model, calibrated with the Kalibr toolbox [17]. Underwater
checkerboard imagery was collected for both cameras to perform the
calibrations. We adapt the camera model code from ORB-SLAM3 [18] to support
our hybrid camera system. In our notation, we represent the camera projection
function as $\pi_{(\cdot)}$, where $\pi_{m}$ and $\pi_{s}$ are the left
monocular stereo camera and rectified stereo projection functions respectively
as defined in [1] and $\pi_{f}$ is the calibrated fisheye projection function
from [16].

###  III-B Feature Representation

While ORB-SLAM2 uses ORB [19] features, ORB performs poorly in many underwater
environments compared to other feature types. We conducted an analysis of the
matching performance of different feature types in the underwater domain,
presented in the results section, which motivated our choice of the SIFT [20]
feature for our system. We adopted CudaSIFT [21], which is currently one of
the fastest GPU accelerated SIFT implementations, for real-time feature
extraction.

###  III-C System Initialization

On system startup, the first keyframe $K_{0}$ is created from the first stereo
frame $F_{s,0}$ that retains at least 8% of the maximum number of features $N$
extracted from the images. $K_{0}$ is set as the origin of the initial map,
which is constructed from all of the stereo keypoints. After the map is
initialized, new keyframes are added from both the stereo and monocular
fisheye frames, with the map scale constrained by the initial and new stereo
map points.

###  III-D Stereo Odometery

Similar to [4], we found that the tracking stage of ORB-SLAM2 failed on our
underwater datasets, even when adapted to SIFT features. A considerable
limitation of the ORB-SLAM2 tracking stage is a constant velocity model, which
has poor accuracy at the low frame rates typical for underwater imaging
systems. [4] used the VISO2 stereo odometer for initial frame pose estimation.
We took inspiration from this and also adopted VISO2 for our system. However,
we found that off-the-shelf VISO2 failed to track our underwater stereo
dataset, due to poor performance of the simple blob and corner response
features, described in the Sobel operator space. We modified VISO2 to use
CudaSIFT features, which are extracted once for each image and then propagated
through the rest of the SLAM pipeline for efficient computation. While the
original VISO2 implementation used a search window to circularly match
features across the current and previous stereo pair, we use GPU accelerated
brute force matching, followed by circular filtering for improved
computational performance. In this scheme, brute force matching is applied
between the previous left and previous right frames, previous right and
current right frames, current right and current left frames, and current left
and previous left frames. A feature is accepted only if the same feature is
matched across all image pairs in a circular fashion. Like in the original
VISO2 implementation, feature matches between a left and right stereo image
pair are further filtered by an epipolar constraint of 1 pixel error
tolerance. However, we found the outlier removal step of the original VISO2 by
2d Delaunay triangulation to be too restrictive in high rugosity coral reef
imagery, resulting in the filtering of many correct feature correspondences.
Through extensive experimentation, we found the circular matching and epipolar
constrained filtering steps were sufficient for removing the majority of
outliers before processing the matches through the ego-motion estimation
stage. The output transform $\widetilde{T}_{s,i}$ of the odometer is the
estimated ego-motion of the current stereo frame $F_{s,i}$ with respect to the
previous stereo frame $F_{s,i-1}$.

###  III-E Tracking

Given an initialized system, the current stereo frame $F_{s,i}$ is processed
through the tracking stage in the following steps:

  1. 1.

Process $F_{s,i}$ through the odometer to extract CudaSIFT features and
estimate the ego-motion $\widetilde{T}_{s,i}$ with respect to the previous
stereo frame.

  2. 2.

Initialize pose of $F_{s,i}$ as $T_{s,i}=\widetilde{T}_{s,i}T_{s,i-1}$, where
$T_{s,i}\in SE(3)$.

  3. 3.

Project the map points tracked in previous frame into the current frame using
projection function $x=\pi_{m}(T_{s,i}X)$, where $X$ is the vector coordinate
of the map point in the world reference frame, and search feature matches
within window regions around the projected points in a brute force manner.

  4. 4.

Formulate the iterative pose optimization:

$\underset{T_{s,i}}{argmin}\sum_{j\in\chi}\rho(\left\|x_{(\cdot)}^{j}-\pi_{(\cdot)}(T_{s,i}X^{j})\right\|_{\Sigma}^{2})$
(1)

where $j\in\chi$ is the set of all map point matches, $\rho$ is the robust
Huber cost function, $\Sigma$ is the covariance matrix associated with the map
point scale, and keypoints can be monocular with $x_{m}^{j}\in R^{2}$ for
projection $\pi_{m}$ or stereo with $x_{s}^{j}\in R^{3}$ for projection
$\pi_{s}$. Only retain inlier map point matches following the optimization.

  5. 5.

If steps 3-4 fail to find enough map point correspondences, track $F_{s,i}$ to
the map points tracked in the reference stereo keyframe $K_{s}$ using BoW
vocabulary levels to guide matching followed by pose optimization as in step
4.

  6. 6.

If few map point matches are found in steps 3-5, set $T_{s,i}$ to the
odometery estimated pose of step 2.

  7. 7.

Search additional map point matches in a local window of keyframes that share
observations of the current map point matches and may include both stereo and
fisheye frames, and repeat the pose optimization of step 4.

If the stereo tracking stage fails completely, the system enters
relocalization mode until tracking is recovered for a stereo frame, as in [1].

In hybrid mode, the current monocular fisheye frame $F_{f,i}$ is processed
through the tracking stage in the following steps, but only after the current
stereo frame is successfully tracked:

  1. 1.

Track $F_{f,i}$ to the map points tracked in the reference fisheye keyframe
$K_{f}$ using BoW vocabulary levels to guide matching followed by pose
optimization as in step 4 of the stereo tracking. The optimization formulation
is the same as Eq. (1), with fisheye projection $\pi_{f}$ and $x_{f}^{j}\in
R^{2}$.

  2. 2.

If step 1 fails, $F_{f,i}$ is tracked to the map points in the reference
stereo keyframe $K_{s}$.

  3. 3.

Search additional map point matches in a local window of keyframes, that may
include both stereo and fisheye frames, and repeat the pose optimization step.

If tracking fails for the fisheye frame but not the current stereo frame, the
system enters relocalization mode for only the fisheye camera, while
continuing mapping of the stereo frames. In this relocalization mode, the
current fisheye frame is first attempted to be matched against all keyframes
in the map using the BoW place recognition to identify match candidates. If
place recognition fails, tracking of the current fisheye frame is then
attempted against the current stereo reference keyframe, proceeding from step
2 above. During the local mapping stage for both stereo and fisheye frames,
the reference keyframe for each is updated to the keyframe that shares the
most feature matches, agnostic to the type of keyframe (i.e. stereo or
fisheye). When a new keyframe is inserted, it is made the reference keyframe
for the next frame of the same type.

During relocalization or when the fisheye frame is tracked against the
reference stereo frame, a perspective-n-point (PnP) solver is constructed to
estimate an initial pose. Like [18], we adopt the Maximum Likelihood
Perspective-n-Point algorithm (MLPnP) [22], which uses projective rays in the
optimization that are agnostic to the camera model, in order to accurately
optimize the feature correspondences between the hybrid fisheye and
perspective stereo frames.

###  III-F Inserting New Keyframes

New stereo and fisheye keyframes are decided following the same scheme as [1]
for stereo and monocular keyframes respectively, with some thresholds tuned
for lower framerates and higher keypoint counts. When a new keyframe is
inserted, new map points are triangulated and added into the map. For each of
these keypoints the maximum and minimum distances that the point can be
detected in a frame are calculated based on the scale of the keypoint in the
reference keyframe. With the hybrid camera system, the scale of the keypoint
can be different at the same distance, depending on which type of frame
observes the keypoint. We resolve this ambiguity by normalizing the keypoint
scale factor by the focal length of the observing frame. This normalization
enables consistent keypoint scale prediction and comparison between hybrid
frames.  As in [1], local bundle adjustment is performed on a set of covisible
keyframes $\mathcal{K}_{L}$ and all points seen in those keyframes
$\mathcal{P}_{L}$ after a new keyframe is inserted into the map, with the
optimization formulated as follows:

$\underset{X^{i},T_{l}}{argmin}\sum_{k\in\mathcal{K}_{L}\cup\mathcal{K}_{F}}\sum_{j\in\chi_{k}}\rho(\left\|x_{(\cdot)}^{j}-\pi_{(\cdot)}(T_{k}X^{j})\right\|_{\Sigma}^{2})$
(2)

where $i\in\mathcal{P}_{L}$, $l\in\mathcal{K}_{L}$, $\chi_{k}$ is the set of
point matches between $\mathcal{P}_{L}$ and keypoints in keyframe $k$, and
$\mathcal{K}_{F}$ is all other keyframes observing $\mathcal{P}_{L}$ but not
in $\mathcal{K}_{L}$ which contribute to the cost but are fixed in the
optimization. Keyframes can be from fisheye images, in which case the
projection function is $\pi_{f}$, or from stereo images in which case the
projection function is $\pi_{m}$ or $\pi_{s}$, depending on whether the
keypoint is monocular or stereo, respectively.

###  III-G Loop Closing

For place recognition, we adapted DBoW2 to SIFT features and trained a million
word vocabulary with ten branching factors and six levels, like the ORB
vocabulary used in ORB-SLAM2. The vocabulary was trained on an extensive set
of underwater imagery data, including the UWHandles and LizardIsland datasets
presented in this paper, plus three large imagery datasets from the Australian
Center for Field Robotics: Tasmania CSP [23], Scott Reef 25 [24], and Tasmania
O’Hara 7 [25]. 2000 CudaSIFT features were extracted per image, with the
images upscaled by a factor of 2 for the first scale pyramid level, the
initial blur set to 1.6, and the difference of Guassian threshold set to 1.0.

After a loop closing event, a full bundle adjustment optimization is
initiated, following the same optimization as Eq. (2), but including all
keyframes and points in the map except for the origin stereo keyframe, which
remains fixed.

###  III-H Datasets

####  III-H1 Stereo Survey Dataset

A stereo SLAM evaluation dataset was collected with a diver operated camera
rig on a shallow coral reef of Lizard Island in Australia (fig. 2). The
dataset was collected using a spiral survey technique [26] that fully covered
a circular area of approximately 14m in diameter, with natural sunlight
providing the only illumination. The rectified stereo image size is 1355x1002
pixels and the images were collected at 5Hz. We refer to this dataset as
LizardIsland.

(a) Left stereo image

(b) COLMAP reconstruction

Figure 2: The LizardIsland spiral survey dataset was collected with a diver
operated stereo rig. The ground truth reconstruction was generated with
COLMAP.

To obtain a ground truth comparison for evaluating our stereo SLAM method, we
processed the dataset through COLMAP [27] to generate a sparse 3D
reconstruction with optimized camera poses. COLMAP does not fix the scale
during optimization, so the reconstruction was scaled in post-process to match
the mean left and right stereo pair baseline to the calibrated value.

####  III-H2 Hybrid Vehicle-Manipulator Dataset

During a cruise in 2019, a hybrid dataset of synchronized vehicle mounted
stereo and wrist mounted fisheye imagery was collected in natural deep ocean
environments of the Costa Rican continental shelf margin with the SuBastian
ROV, operated by Schmidt Ocean Institute [28]. We have previously published
the fisheye imagery portion of this data as the UWHandles dataset [29]. For
this work, we have extended this dataset by further processing four
environmentally unique stereo and fisheye image sequences for evaluation of
our hybrid SLAM method. We refer to these sequences as Mounds1, Mounds2,
Seeps1, and Seeps2. For these sequences, TagSLAM [30] was used to obtain
ground truth pose estimates for the stereo and fisheye cameras, based on the
detection of AprilTags [31] distributed in the scenes.

##  IV Evaluation

All evaluations were run on a desktop computer with an AMD Ryzen Threadripper
2990WX CPU and an NVIDIA Titan V GPU.

###  IV-A Comparative Feature Analysis

We conducted an evaluation to determine which feature representation is best
adapted to the visual degradation of underwater environments and can be
robustly matched between hybrid perspective and fisheye frames with variable
relative poses. We sampled every fifth hybrid frame from each of the UWHandles
dataset sequences and, to reduce any bias from artificial features, we used
the tracked AprilTag poses from TagSLAM to project circular masks over the
tags in the image frames. For each feature type, 2000 features were extracted
from each image, and the features were brute force matched across each hybrid
fisheye and left stereo image pair. Lowe’s ratio test was applied to remove
ambiguous matches, with a ratio threshold of 0.8 for all feature types except
ContextDesc, which achieved significantly improved performance with a ratio of
0.9. OpenCV’s RANSAC based essential matrix fitting was used to filter the
matches and recover a relative pose estimate between each fisheye and stereo
frame. Table I shows the results of this evaluation. Given that an essential
matrix based pose estimate does not provide scale, both the orientation and
translation errors of the pose estimates were evaluated as angular errors. For
translation, this error is the angular difference between the translation
direction vector from the left stereo frame to the fisheye frame. The
performance was evaluated using the area under the accuracy-threshold curve
(AUC) with a max angular error of 180°. While most of the tested feature types
were popular conventional features, we also tested two deep learned feature
variants: ContextDesc and SuperPoint. We note that the learned features were
used with their provided model weights and were not fine-tuned on underwater
data. Of the conventional feature types, ROOT_SIFT and SIFT perform the best,
achieving significantly better performance than ORB. Of the deep learned
features, SuperPoint had highly variable performance across the different
sequences, and the mean number of inlier matches was lower than other
conventional features. Interestingly, ContextDesc performed the best overall
out of all the feature types, consistently matching more than double the
features of ROOT_SIFT and achieving very high AUC scores. It is noteworthy
that ContextDesc uses SIFT interest points but learns the descriptor, so all
of the best performing features are based on the SIFT detector. These results
merit further investigation into the application of learned features for
underwater vision. For our initial implementation in this work, we chose to
use a highly optimized GPU accelerated implementation of SIFT, but we note
that the learned descriptors of ContextDesc are 128-d, like SIFT, and are
directly compatible with our entire method pipeline.

TABLE I: Area under accuracy-threshold curve evaluation of feature matching performance on the UWHandles underwater hybrid image sequences. Accuracy is evaluated as angular error in the predicted rotation (AUC Rot) and translation direction vector (AUC Trans) between each hybrid left stereo and fisheye image pair. Also reported is the mean number of inlier feature matches across each sequence. Sequence |  |  SIFT [20] |  ROOT SIFT[32] |  ORB [33] |  SURF [34] |  AKAZE [35] |  CONTEXTDESC [36] |  SUPERPOINT [37]
---|---|---|---|---|---|---|---|---
Mounds1 | AUC Trans | 0.949 | 0.936 | 0.856 | 0.877 | 0.922 | 0.970 | 0.98
AUC Rot | 0.937 | 0.948 | 0.750 | 0.812 | 0.858 | 0.951 | 0.964
Mean Matches | 91 | 101 | 25 | 34 | 46 | 210 | 74
Mounds2 | AUC Trans | 0.946 | 0.940 | 0.864 | 0.886 | 0.906 | 0.976 | 0.947
AUC Rot | 0.810 | 0.853 | 0.488 | 0.676 | 0.629 | 0.959 | 0.873
Mean Matches | 31 | 33 | 14 | 21 | 21 | 80 | 41
Seeps1 | AUC Trans | 0.964 | 0.980 | 0.885 | 0.869 | 0.904 | 0.986 | 0.938
AUC Rot | 0.944 | 0.965 | 0.745 | 0.770 | 0.811 | 0.980 | 0.885
Mean Matches | 64 | 73 | 23 | 28 | 43 | 150 | 53
Seeps2 | AUC Trans | 0.960 | 0.964 | 0.935 | 0.930 | 0.954 | 0.974 | 0.917
AUC Rot | 0.942 | 0.953 | 0.894 | 0.891 | 0.926 | 0.965 | 0.763
Mean Matches | 89 | 100 | 60 | 50 | 92 | 146 | 39

(a) 4k features/frame

(b) 2k features/frame

Figure 3: Final stereo SLAM maps on the LizardIsland dataset, showing the
densely connected keyframe graphs.

###  IV-B Stereo SLAM

(a) Tracking

(b) SLAM

Figure 4: Stereo SLAM results (green) for the LizardIsland dataset when loop
closing is disabled (4(a)) and enabled (4(b)), plotted against the ground
truth (black).

The core of our system is a stereo SLAM pipeline, which must be robust to
underwater environments. We used the LizardIsland survey dataset to evaluate
the stereo SLAM performance. We tested ORB-SLAM2 on this dataset, both with
and without loop closing enabled, but it lost track after only a few frames
and was unable to relocalize. We also tested the vanilla VISO2 stereo
odometer, but, even with extensive tuning, VISO2 failed to track the dataset
with sensible accuracy. As a result, the stereo SLAM method of [4],
specifically designed for underwater environments, also failed on our dataset,
as it runs on top of VISO2. We evaluated our stereo SLAM method with both 2000
and 4000 CudaSIFT features extracted each frame, with an interest point
Difference of Gaussian threshold of 1.2. Figure 4 shows the results for 4000
features, both with and without loop closing enabled, and table II gives the
performance of the system in all tests. The test trajectories were aligned
with the COLMAP ground truth using the Horn method [38] without scaling. The
figure shows that our visual odometer based tracking method without loop
closing tracks very well in the horizontal plane with most of the drift error
being accumulated in the z-depth estimate. For both extracted feature counts,
the table shows that a high accuracy, with less than 2cm root mean squared
absolute trajectory error (RMSE), is attained by the full SLAM system with
loop closing. For 4000 features per frame, the number of map points in the
final map is approximately double the map point count for 2000 features per
frame, showing that the system scales well with the number of extracted
features. The system can achieve >10Hz for 2000 features per frame, which is a
high framerate for underwater systems. Figure 3 shows the densely connected
keyframe graphs for the final SLAM maps, demonstrating consistent loop closing
between neighboring spiral trajectories.

TABLE II: Stereo SLAM performance on the LizardIsland dataset, with the number of extracted features is set to 4000 and 2000. Performance is evaluated as RMSE of the absolute trajectory error. Results are reported with and without loop closing enabled. Also reported is the number of keyframes (KFs) and map points (MPs) in the final map and the average frame processing time in the tracking thread. System Mode | RMSE (cm) | KFs | MPs | Avg Time (ms)
---|---|---|---|---
Tracking Only (4000) | 49.1 | - | - | 94.2
Loop Closing (4000) | 1.4 | 562 | 190,474 | 117.9
Tracking Only (2000) | 58.2 | - | - | 55.7
Loop Closing (2000) | 1.8 | 622 | 98,812 | 64.8

###  IV-C Hybrid SLAM

TABLE III: Hybrid SLAM timing evaluation, measured as the mean frame processing time in the tracking thread. # Features / Frame | 2k | 4k | 6k
---|---|---|---
Mean Time | 179ms | 249ms | 314ms

We evaluated the performance of our hybrid SLAM system on the four sequences
of the UWHandles dataset. The results are reported in table IV. For all
sequences, every stereo frame was successfully registered in the SLAM map.
Given that the stereo camera is mostly stationary across these image
sequences, and to reduce the effect of noise in the imperfect ground truth, we
evaluated the hybrid SLAM error using the relative pose estimates between the
left stereo and fisheye cameras for each synchronized hybrid frame. Only
hybrid frames where the fisheye frame was successfully registered into the map
were included in the error evaluation. As is detailed in the table, the system
generated approximately twice as many keyframes and map points when running in
hybrid mode versus stereo only mode, demonstrating the ability to extend the
map beyond the limited stereo camera viewpoint. Also, the majority of fisheye
frames were successfully registered into the map for all sequences. Despite
the sequences varying significantly in environment type, the hybrid SLAM mode
is able to generate a similar amount of keyframes and map points for each
sequence, and the estimated pose errors are very similar across each sequence,
demonstrating our system can operate in challenging and diverse, natural
seafloor environments. Figure 5 shows a frame capture from running hybrid SLAM
on each of the sequences. Table III provides the timing evaluation for
processing a hybrid stereo and fisheye frame pair through the tracking thread
for different feature count settings. For 4000 features extracted per image,
the system can easily attain 3Hz, which is the rate that the UWHandles data
was collected.

(a) Mounds1

(b) Mounds2

(c) Seeps1

(d) Seeps2

Figure 5: Snapshots of hybrid SLAM running on the UWHandles sequences. Top row is the left stereo camera frame, middle row is the manipulator mounted fisheye frame, and bottom row is the map with the keypoints and keyframes. TABLE IV: Evaluation of hybrid SLAM on the UWHandles dataset. Error is evaluated on the estimated pose difference between the left stereo and fisheye cameras for each synchronized hybrid frame, where $\Delta t$ is translation error and $\Delta q$ is rotation error. The ”hybrid matches” column gives the number of fisheye frames registered in the map over the total number of frames in the sequence. The error is only evaluated over the registered frames. The ”KFs” column is the number of keyframes in the final map for the hybrid SLAM mode versus stereo only mode, and the ”MPs” column is the same format for the number of final keypoints in the map. Sequence |  $\Delta t$ mean (cm) |  $\Delta t$ median (cm) |  $\Delta q$ mean (deg) |  $\Delta q$ median (deg) | hybrid matches | KFs hybrid/stereo | MPs hybrid/stereo
---|---|---|---|---|---|---|---
Mounds1 | 2.04 | 2.06 | 0.58 | 0.50 | 652 / 783 | 21 / 11 | 4271 / 2671
Mounds2 | 1.38 | 1.22 | 0.98 | 0.82 | 713 / 756 | 24 / 11 | 4086 / 1773
Seeps1 | 2.82 | 2.02 | 1.17 | 0.46 | 1059 / 1089 | 24 / 11 | 4636 / 2847
Seeps2 | 2.12 | 1.84 | 1.38 | 0.97 | 778 / 802 | 23 / 16 | 4320 / 2365

##  V Conclusion

We have presented a novel hybrid SLAM method, targeting deployment on
underwater vehicle manipulator systems, that can operate in real-time. The
method can fuse features from both a vehicle mounted stereo camera and a
manipulator mounted fisheye camera into the same map, enabling dynamic
viewpoint acquisition and map extension with the manipulator mounted camera.
We have demonstrated the robustness of our method on both a shallow reef
stereo image survey dataset and on four hybrid image sequences captured in
natural, deep seafloor environments. In this work, the system has been tested
on a desktop computer, which is suitable for ROV operations, where the compute
can run topside through a tethered connection with the vehicle. In future
work, we will optimize the method to run on embedded GPU devices which can be
deployed onboard the vehicle. Currently, the system fuses the fisheye data
into the map without consideration of the manipulator state. In future work,
we will also explore the formulation of a kinematic factor from the
manipulator joint states between the fisheye and the stereo camera to improve
registration of the fisheye camera into the map and the overall robustness of
the mapping method. This factor would also enable real-time feedback for the
kinematic calibration of the manipulator, which is a challenging problem for
the imprecise hydraulic manipulators common for underwater systems. We will
also explore the use of learned feature descriptors such as ContextDesc to
improve system performance. Finally, we will explore methods for generating
dense reconstructions based on the sparse feature maps and optimized camera
poses to build a complete real-time scene reconstruction method for UVMSs.

## References

  * [1] Raul Mur-Artal and Juan D Tardós  “Orb-slam2: An open-source slam system for monocular, stereo, and rgb-d cameras”  In _IEEE Transactions on Robotics_ 33.5 IEEE, 2017, pp. 1255–1262
  * [2] Danping Zou and Ping Tan  “Coslam: Collaborative visual slam in dynamic environments”  In _IEEE transactions on pattern analysis and machine intelligence_ 35.2 IEEE, 2012, pp. 354–366
  * [3] M. Johnson-Roberson et al.  “High-resolution Underwater Robotic Vision-based Mapping and 3D Reconstruction for Archaeology”  In _Journal of Field Robotics_ , 2016
  * [4] Pep Lluis Negre, Francisco Bonin-Font and Gabriel Oliver  “Cluster-based loop closing detection for underwater slam in feature-poor regions”  In _2016 IEEE International Conference on Robotics and Automation (ICRA)_ , 2016, pp. 2589–2595  DOI: 10.1109/ICRA.2016.7487416
  * [5] Maxime Ferrera, Julien Moras, Pauline Trouvé-Peloux and Vincent Creuze  “Real-time Monocular Visual Odometry for Turbid and Dynamic Underwater Environments” arXiv: 1806.05842  In _arXiv:1806.05842 [cs]_ , 2020  URL: http://arxiv.org/abs/1806.05842
  * [6] Franco Hidalgo, Chris Kahlefendt and Thomas Bräunl  “Monocular ORB-SLAM Application in Underwater Scenarios”  In _2018 OCEANS-MTS/IEEE Kobe Techno-Oceans (OTO)_ , 2018, pp. 1–4  IEEE
  * [7] Paul Ozog, Matthew Johnson-Roberson and Ryan M Eustice  “Mapping underwater ship hulls using a model-assisted bundle adjustment framework”  In _Robotics and Autonomous Systems_ 87 Elsevier, 2017, pp. 329–347
  * [8] Andreas Geiger, Julius Ziegler and Christoph Stiller  “StereoScan: Dense 3D Reconstruction in Real-time”  In _Intelligent Vehicles Symposium (IV)_ , 2011
  * [9] Pep Lluis Negre Carrasco, Francisco Bonin-Font and Gabriel Oliver-Codina  “Global image signature for visual loop-closure detection”  In _Autonomous Robots_ 40.8 Springer, 2016, pp. 1403–1417
  * [10] Matthew Klingensmith, Siddartha S Sirinivasa and Michael Kaess  “Articulated robot motion for simultaneous localization and mapping (arm-slam)”  In _IEEE robotics and automation letters_ 1.2 IEEE, 2016, pp. 1156–1163
  * [11] Jinghui Li, Akitoshi Ito and Yusuke Maeda  “A SLAM-Integrated Kinematic Calibration Method for Industrial Manipulators with RGB-D Cameras”  In _2019 19th International Conference on Control, Automation and Systems (ICCAS)_ , 2019, pp. 686–689  IEEE
  * [12] Mingyou Chen et al.  “3D global mapping of large-scale unstructured orchard integrating eye-in-hand stereo vision and SLAM”  In _Computers and Electronics in Agriculture_ 187 Elsevier, 2021, pp. 106237
  * [13] Arun Das and Steven L. Waslander  “Calibration of a dynamic camera cluster for multi-camera visual SLAM”  In _2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2016, pp. 4637–4642  DOI: 10.1109/IROS.2016.7759682
  * [14] Dorian Gálvez-López and J. D. Tardós  “Bags of Binary Words for Fast Place Recognition in Image Sequences”  In _IEEE Transactions on Robotics_ 28.5, 2012, pp. 1188–1197  DOI: 10.1109/TRO.2012.2197158
  * [15] Morgan Quigley et al.  “ROS: an open-source Robot Operating System”  In _ICRA workshop on open source software_ 3.3.2, 2009, pp. 5  Kobe, Japan
  * [16] Juho Kannala and Sami S Brandt  “A generic camera model and calibration method for conventional, wide-angle, and fish-eye lenses”  In _IEEE transactions on pattern analysis and machine intelligence_ 28.8 IEEE, 2006, pp. 1335–1340
  * [17] Jérôme Maye, Paul Furgale and Roland Siegwart  “Self-supervised calibration for robotic systems”  In _2013 IEEE Intelligent Vehicles Symposium (IV)_ , 2013, pp. 473–480  IEEE
  * [18] Carlos Campos et al.  “ORB-SLAM3: An Accurate Open-Source Library for Visual, Visual–Inertial, and Multimap SLAM”  In _IEEE Transactions on Robotics_ IEEE, 2021
  * [19] Ethan Rublee, Vincent Rabaud, Kurt Konolige and Gary Bradski  “ORB: An efficient alternative to SIFT or SURF”  In _2011 International conference on computer vision_ , 2011, pp. 2564–2571  Ieee
  * [20] David G Lowe  “Distinctive image features from scale-invariant keypoints”  In _International journal of computer vision_ 60.2 Springer, 2004, pp. 91–110
  * [21] Mårten Björkman, Niklas Bergström and Danica Kragic  “Detecting, segmenting and tracking unknown objects using multi-label MRF inference”  In _Computer Vision and Image Understanding_ 118 Elsevier, 2014, pp. 111–127
  * [22] Steffen Urban, Jens Leitloff and S. Hinz  “MLPNP – A REAL-TIME MAXIMUM LIKELIHOOD SOLUTION TO THE PERSPECTIVE-N-POINT PROBLEM”  In _ISPRS Annals of Photogrammetry, Remote Sensing and Spatial Information Sciences_ III-3, 2016, pp. 131–138  DOI: 10.5194/isprsannals-III-3-131-2016
  * [23] Lisa Meyer, Nicole Hill, Peter Walsh and Neville Barrett  “Methods for the processing and scoring of AUV digital imagery from South Eastern Tasmania”  In _Institute for Marine and Antarctic Studies Internal Report_ Citeseer, 2011
  * [24] Daniel M Steinberg, Stefan B Williams, Oscar Pizarro and Michael V Jakuba  “Towards autonomous habitat classification using Gaussian Mixture Models”  In _2010 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , 2010, pp. 4424–4431  IEEE
  * [25] D Steinberg, Ariell Friedman, Oscar Pizarro and Stefan B Williams  “A Bayesian nonparametric approach to clustering data from underwater robotic surveys”  In _International Symposium on Robotics Research_ 28, 2011, pp. 1–16  Citeseer
  * [26] Oscar Pizarro et al.  “A simple, fast, and repeatable survey method for underwater visual 3D benthic mapping and monitoring”  In _Ecology and Evolution_ 7.6 Wiley Online Library, 2017, pp. 1770–1782
  * [27] Johannes Lutz Schönberger and Jan-Michael Frahm  “Structure-from-Motion Revisited”  In _Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2016
  * [28] Peter Vrolijk et al.  “Using a Ladder of Seeps With Computer Decision Processes to Explore for and Evaluate Cold Seeps on the Costa Rica Active Margin”  In _Frontiers in Earth Science_ 9 Frontiers, 2021, pp. 143
  * [29] Gideon Billings and Matthew Johnson-Roberson  “SilhoNet-fisheye: Adaptation of a ROI based object pose estimation network to monocular fisheye images”  In _IEEE Robotics and Automation Letters_ 5.3 IEEE, 2020, pp. 4241–4248
  * [30] Bernd Pfrommer and Kostas Daniilidis  “Tagslam: Robust slam with fiducial markers”  In _arXiv preprint arXiv:1910.00679_ , 2019
  * [31] John Wang and Edwin Olson  “AprilTag 2: Efficient and robust fiducial detection”  In _2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2016, pp. 4193–4198  IEEE
  * [32] Relja Arandjelović and Andrew Zisserman  “Three things everyone should know to improve object retrieval”  In _2012 IEEE Conference on Computer Vision and Pattern Recognition_ , 2012, pp. 2911–2918  IEEE
  * [33] Ethan Rublee, Vincent Rabaud, Kurt Konolige and Gary Bradski  “ORB: An efficient alternative to SIFT or SURF” ISSN: 2380-7504  In _2011 International Conference on Computer Vision_ , 2011, pp. 2564–2571  DOI: 10.1109/ICCV.2011.6126544
  * [34] Herbert Bay, Andreas Ess, Tinne Tuytelaars and Luc Van Gool  “Speeded-up robust features (SURF)”  In _Computer vision and image understanding_ 110.3 Elsevier, 2008, pp. 346–359
  * [35] Pablo F Alcantarilla and T Solutions  “Fast explicit diffusion for accelerated features in nonlinear scale spaces”  In _IEEE Trans. Patt. Anal. Mach. Intell_ 34.7, 2011, pp. 1281–1298
  * [36] Zixin Luo et al.  “Contextdesc: Local descriptor augmentation with cross-modality context”  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 2527–2536
  * [37] Daniel DeTone, Tomasz Malisiewicz and Andrew Rabinovich  “Superpoint: Self-supervised interest point detection and description”  In _Proceedings of the IEEE conference on computer vision and pattern recognition workshops_ , 2018, pp. 224–236
  * [38] Berthold KP Horn  “Closed-form solution of absolute orientation using unit quaternions”  In _Josa a_ 4.4 Optical Society of America, 1987, pp. 629–642


# Coherent Electric-Field Control of Orbital State in a Neutral Nitrogen-
Vacancy Center

Hodaka Kurokawa<EMAIL_ADDRESS>Quantum Information Research
Center, Institute of Advanced Sciences, Yokohama National University, 79-5
Tokiwadai, Hodogaya, Yokohama 240-8501, Japan  Keidai Wakamatsu  Department of
Physics, Graduate School of Engineering Science, Yokohama National University,
79-5 Tokiwadai, Hodogaya, Yokohama 240-8501, Japan  Shintaro Nakazato
Department of Physics, Graduate School of Engineering Science, Yokohama
National University, 79-5 Tokiwadai, Hodogaya, Yokohama 240-8501, Japan
Toshiharu Makino  Quantum Information Research Center, Institute of Advanced
Sciences, Yokohama National University, 79-5 Tokiwadai, Hodogaya, Yokohama
240-8501, Japan  Advanced Power Electronics Research Center, National
Institute of Advanced Industrial Science and Technology, 1-1-1 Umezono,
Tsukuba, Ibaraki, 305-8568, Japan  Hiromitsu Kato  Quantum Information
Research Center, Institute of Advanced Sciences, Yokohama National University,
79-5 Tokiwadai, Hodogaya, Yokohama 240-8501, Japan  Advanced Power Electronics
Research Center, National Institute of Advanced Industrial Science and
Technology, 1-1-1 Umezono, Tsukuba, Ibaraki, 305-8568, Japan  Yuhei Sekiguchi
Quantum Information Research Center, Institute of Advanced Sciences, Yokohama
National University, 79-5 Tokiwadai, Hodogaya, Yokohama 240-8501, Japan  Hideo
Kosaka<EMAIL_ADDRESS>Quantum Information Research Center,
Institute of Advanced Sciences, Yokohama National University, 79-5 Tokiwadai,
Hodogaya, Yokohama 240-8501, Japan  Department of Physics, Graduate School of
Engineering Science, Yokohama National University, 79-5 Tokiwadai, Hodogaya,
Yokohama 240-8501, Japan

###### Abstract

The coherent control of the orbital state is crucial for color centers in
diamonds for realizing extremely low-power manipulation. Here, we propose the
neutrally charged nitrogen-vacancy center, NV0, as an ideal system for orbital
control through electric fields. We estimate electric susceptibility in the
ground state of NV0 to be comparable to that in the excited state of NV-.
Also, we demonstrate coherent control of the orbital states of NV0. The
required power for orbital control is three orders of magnitude smaller than
that for spin control, highlighting the potential for interfacing a
superconducting qubit operated in a dilution refrigerator.

Diamond color centers are attracting broad attention because of their
potential applications in quantum communication [1, 2, 3], quantum computation
[4, 5], and quantum sensing [6, 7]. The spin degree of freedom is primarily
utilized as a quantum bit, thanks to its long coherence time of over one
second [8, 9, 10], and excellent controllability [11, 12]. On the other hand,
control over the orbital degree of freedom is also crucial for some
applications such as frequency tuning of the zero-phonon line photons and
extremely low-power control of the electron states. The ability to tune the
zero-phonon line frequency through electric fields or strain is essential for
generating entanglement between remote color centers [13, 14, 1]. Moreover,
the coupling of electric fields or strain with the orbital degree of freedom
is more stronger compared to the magnetic field’s coupling with the spin [15,
16, 17], enabling highly efficient control of the electron state. With strong
spin-orbit coupling, efficient spin-state control via strain has been achieved
in color centers [18], which is particularly advantageous for operation in a
dilution refrigerator. Nevertheless, directly achieving coherent control over
the orbital state remains challenging for representative color centers due to
the short lifetime of the optically excited state of NV- ($\sim 10$ ns) [19]
and large ground-state splitting of group-IV color centers [20].

Here, we propose the neutrally charged nitrogen-vacancy center, NV0 [21, 22,
23, 24, 25, 26, 27, 28], as an ideal system for orbital state control through
electric fields. The ground-state spin-orbit splitting of NV0 is approximately
10 GHz [25, 26], enabling direct microwave access. Additionally, the ground
state of NV0 exhibits an orbital relaxation time of several hundreds of
nanoseconds [26], which is more than one order of magnitude longer than that
of the excited state of NV-. Furthermore, the energy level structure of the
ground state in NV0 is similar to that of group-IV color centers. Therefore,
NV0 can serve as an ideal system for understanding the properties of the
orbital state and conducting proof-of-principle experiments relevant to group-
IV color centers. In this paper, we investigate the electric susceptibility of
NV0 and demonstrate coherent control over its orbital state. The highly
efficient control opens up possibilities for future applications, such as a
quantum interface communicating with a superconducting quantum bit in a
dilution refrigerator [29, 30].

## Results

### Characterization of NV0

Figure 1: The energy level structure of NV0, the device image, the
photoluminescence excitation (PLE) spectrum, and the dc voltage dependence of
the PLE spectrum. (a) Schematic of the energy level of NV0 with some static
strain under a zero magnetic field. ground state exhibits an energy splitting
on the order of 10 GHz due to the combined effects of spin-orbit interaction
and strain. The optical transition occurs at 575 nm. (b) An optical microscope
image of the electrical circuit on the diamond used for experiments. The upper
electrodes are used to apply low ($<$ 1 MHz) or high ($\sim$10 GHz) electric
fields. The color map shows the photoluminescence (PL) counts around the
center of the circuit. The red-circled NV center is used for experiments. (c)
The PLE spectrum of NV0 and its measurement sequence. The 637 nm laser
resonant to the zero-phonon line of NV- is used to initialize the charge state
to NV0. The 575 nm laser is used to observe the transition. (d) The PLE
spectrum shifts as a function of dc voltage. Green and blue circles
corresponds to the lower ($\ket{0}$) and upper ($\ket{1}$) branches,
respectively. The inset shows the difference of the PLE frequencies between
upper ($\ket{0}$) and lower $\ket{1}$ branches. The horizontal axis is the
same as the main graph.

Under a zero-magnetic field, the ground-state Hamiltonian, $H$, of NV0 in
$\ket{\pm}_{\text{o}}$ basis can be modeled as [25, 26] :

$H^{(\pm)}/h=2\lambda\hat{L}_{z}\hat{S}_{z}+\epsilon_{\perp}(\hat{L}_{+}+\hat{L}_{-}),$
(1)

where $h$ is the Plank’s constant, $\lambda$ is the spin-orbit interaction
parameter, $\hat{L}_{z}=\sigma_{z}$,
$\hat{L}_{\pm}=\ket{\pm}_{\text{o}}\bra{\mp}_{\text{o}}$ are the orbital
operator in $\ket{\pm}_{\text{o}}=\mp 1/\sqrt{2}(\ket{e_{x}}_{\text{o}}\pm
i\ket{e_{y}}_{\text{o}})$ basis, $\ket{e_{x}}_{\text{o}}$ and
$\ket{e_{y}}_{\text{o}}$ are the strain eigenstates, $S_{z}=(1/2)\sigma_{z}$
is the spin operator for the 1/2 spin, and $\epsilon_{\perp}$ is the
perpendicular strain parameter. Here, the subscripts o (s) denote the orbital
(spin) degree of freedom. Figure 1 (a) shows the energy level of NV0. The
ground-state splitting comes from the spin-orbit coupling, $\lambda$, and the
strain, $\epsilon_{\perp}$, as $2\sqrt{\lambda^{2}+\epsilon_{\perp}^{2}}$.
Under the strain, the eigenstates change from $\ket{\pm}_{\text{o}}$ to
$\ket{+^{\prime}}_{\text{o}}=\alpha\ket{+}_{\text{o}}+\beta\ket{-}_{\text{o}}$
and
$\ket{-^{\prime}}_{\text{o}}=-\beta\ket{+}_{\text{o}}+\alpha\ket{-}_{\text{o}}$
with $\alpha^{2}+\beta^{2}=1$. Figure 1 (b) shows the electrical circuit
directly formed on a diamond. The low ($<$ 1 MHz) and high ($\sim$10 GHz)
frequency electric fields are applied to nitrogen-vacancies using the upper
electrodes. The color map of photoluminescence (PL) counts around the center
of the circuit is also shown. Figure 1 (c) shows the photoluminescence
excitation (PLE) spectrum with the experimental sequence. A 637 nm red laser
resonant to zero-phonon line (ZPL) of NV- is used to convert NV- to NV0.
Frequencies of a 575 nm yellow laser are swept to search ZPL of NV0. Two
transition lines are observed, indicating that the ground-state splitting is
12.85 GHz.

Since all the experiments are performed under an ambient magnetic field, the
transitions from $\ket{+^{\prime}}_{\mathrm{o}}\ket{\uparrow}_{\mathrm{s}}$
and $\ket{-^{\prime}}_{\mathrm{o}}\ket{\downarrow}_{\mathrm{s}}$ can not be
distinguished. Also,
$\ket{-^{\prime}}_{\mathrm{o}}\ket{\uparrow}_{\mathrm{s}}$ and
$\ket{+^{\prime}}_{\mathrm{o}}\ket{\downarrow}_{\mathrm{s}}$ are
indistinguishable. In the following, we write as
$\ket{-^{\prime}}_{\mathrm{o}}\ket{\uparrow}_{\mathrm{s}}=\ket{+^{\prime}}_{\mathrm{o}}\ket{\downarrow}_{\mathrm{s}}\equiv\ket{0}$,
$\ket{+^{\prime}}_{\mathrm{o}}\ket{\uparrow}_{\mathrm{s}}=\ket{-^{\prime}}_{\mathrm{o}}\ket{\downarrow}_{\mathrm{s}}\equiv\ket{1}$,
and
$\ket{0}_{\mathrm{o}}\ket{\uparrow}_{\mathrm{s}}=\ket{0}_{\mathrm{o}}\ket{\downarrow}_{\mathrm{s}}\equiv\ket{2}$
ignoring the spin degree of freedom. In other words, we treat the four-level
system in the ground state as effective two-level systems with different spin
states. Since the lifetime of the spin reaches several milliseconds even under
laser irradiation [26], the assumption is held as long as the laser
irradiation time is sufficiently short.

Next, we apply dc electric fields to estimate the parameters in eq. (1) and
the electric susceptibility. The effects of both strain and dc electric fields
in $\ket{\pm}_{\text{o}}\ket{\uparrow}_{\text{s}}$ basis can be written as:

$\begin{split}H^{(\pm)}/h=\lambda\hat{L}_{z}+d_{\parallel}E_{z}\hat{I}+(\epsilon_{\perp}+d_{\perp}E_{\perp})(\hat{L}_{+}+\hat{L}_{-})\\\
+d_{\perp}E_{\perp}^{\prime}(-i\hat{L}_{+}+i\hat{L}_{-}),\end{split}$ (2)

where $d_{\parallel}$ ($d_{\perp}$) is the electric susceptibility parallel
(perpendicular) to the NV axis, $E_{z}$ is the electric fields parallel to the
NV axis, $E_{\perp}$ is the electric fields whose axis is same with the
direction of $\epsilon_{\perp}$, $E_{\perp}^{\prime}$ is the electric fields
perpendicular to $E_{\perp}$ and $E_{z}$. Here, we ignore the strain term
which only contributes to the global shift of the energy (See Supplementary
materials for details). The energy eigenvalues are
$d_{\parallel}E_{z}\pm\sqrt{\lambda^{2}+(\epsilon_{\perp}+E_{\perp})^{2}+(d_{\perp}E_{\perp}^{\prime})^{2}}$.
Figure 1 (d) shows the PLE frequency shifts of both upper ($\ket{1}$) and
lower ($\ket{0}$) branches as a function of the dc electric fields. With
increasing the dc electric fields, the PLE frequencies of the two branches
decrease almost linearly. With the application of $\pm$50 V, the PLE
frequencies can be shifted up to $\sim$30 GHz. The linear frequency shifts in
both branches are caused by $E_{z}$. Using the electric fields distribution
obtained from the finite element simulation, $d_{\parallel}$ is estimated to
be 1.08 MHz/(V cm-1) (See Supplementary materials for details), which is
similar to the reported value in the excited state of NV-, 0.7 MHz/(V cm-1)
[17]. Additionally, the frequency shifts due to $E_{\perp}$ and
$E_{\perp}^{\prime}$ can be seen clearly by subtracting the upper PLE
frequencies from the lower PLE frequencies (The inset of Fig. 1 (d)). The
frequency differences are fitted with the equation,
$2\sqrt{\lambda^{2}+(\epsilon_{\perp}+d_{\perp}E_{\perp})^{2}+(d_{\perp}E_{\perp}^{\prime})^{2}}$
with a fixed $\lambda=4.80$ GHz [26]. The parameters are estimated to be
$\epsilon_{\perp}$ = 4.06 GHz and $d_{\perp}$ = 363 kHz/(V cm-1). Note that
$d_{\perp}$ is several times smaller than the value reported in Ref. [17] for
the excited state of NV-, 1.4 MHz/(V cm-1).

Figure 2: Measurement of the orbital relaxation time of NV0. (a) Schematic of
the experimental sequence for the measurement of the orbital relaxation time
of NV0. After the charge initialization (Init.) using the 637 nm laser, the
first 575 nm laser pulse initializes $\ket{0}$ to $\ket{1}$ and the population
in $\ket{0}$ is read out (RO) by the second 575 nm pulse after some delay
time. The peak counts, $A$ and $B$, at each pulse are used to estimate the
population in $\ket{0}$. The repetition times for each measurement, $N$, are
150 which is determined from the measurement of the decay of the PL counts as
a function of the duration of the 575 nm laser irradiation. (b) The normalized
counts, $B/A$, as a function of the delay time. The red dots are the
experimental data and the orange curve is the result of the curve fit with the
exponential function. The inset shows The time-resolved PL counts with 100 ns
delay. The PL counts are smoothed using the Savitzky-Golay filter with 17
points before the calculation of the peak height. The background is subtracted
using the data when the laser is absent.

### The orbital relaxation of NV0

The orbital relaxation time, $T_{1}$, is important since it determines the
limit of operation time. Figure 2 (a) shows the experimental sequence for the
measurement of $T_{1}$. A 1-$\mu$s long optical pulse resonant to the
transition between $\ket{0}$ and $\ket{2}$ initializes NV0 into $\ket{1}$ from
thermal mixture of $\ket{0}$ and $\ket{1}$ (The inset of Fig. 2 (b)). The
first peak height is proportional to the initial population in $\ket{0}$. The
peak height in the second pulse is proportional to the population in $\ket{0}$
relaxed from $\ket{1}$. The ratio of the pulse heights in the first and second
pulses is a measure of decay from $\ket{1}$, which are used to calculate
$T_{1}$. Figure 2 (b) shows the normalized pulse height as a function of the
delay of the second pulse. $T_{1}$ is estimated to be $\sim 137$ ns at 5.5 K
from the curve fit with $1-a\mathrm{exp}(-t/T_{1})+b$, where $a$ and $b$ are
some constants. In Ref. [26], $T_{1}$ is a few times larger than ours. At 5.5
K, $T_{1}$ should be limited by thermal phonons. Further increase of $T_{1}$
is expected by decreasing temperature down to several tens of millikelvin.

Figure 3: The optically detected electrical resonance (ODER) and the Rabi
splitting measurements. (a) Schematic of the experimental sequence for the
measurements of ODER and the Rabi splitting. After the charge conversion by
the 637 nm laser, the 575 nm laser resonant to the transition between
$\ket{0}$ and $\ket{2}$ is applied simultaneously with microwave application.
Frequencies of the 575 nm laser are swept for the ODER measurement. The power
of the microwave is swept for the measurement of the Rabi splitting. (b) The
ODER spectrum of NV0. The orange dots are the data and the red curve is the
Gaussian fit. (c) The Rabi splitting as a function of the square root of the
power of the applied microwave electric fields.

### Optically detected electrical resonance and the Rabi splitting

To investigate the electrical resonance frequency in the ground state of NV0,
we employ the optically detected electrical resonance (ODER). Since we apply
the driving electric fields resonant to the eigenfrequency determined by
$\lambda$ and $\epsilon_{\perp}$, we treat the ac driving electric fields as a
perturbation to the eigenstates of the strain,
$\ket{\pm^{\prime}}_{\text{o}}$. In the rotating frame of the driving ac
electric fields, $H_{\mathrm{d}}/h=f_{\mathrm{d}}\hat{L}_{z}/2$, under the
rotating wave approximation, Hamiltonian of NV0 in
$\ket{\pm^{\prime}}_{\textrm{o}}$ basis can be written as:

$H^{(\pm^{\prime})}/h=\Delta\hat{L}_{z}+d_{\perp}E_{\perp}^{\prime\prime}(\hat{L}_{+}+\hat{L}_{-}),$
(3)

where $f_{\text{d}}$ is the frequency of the driving electric fields,
$\Delta=\sqrt{\lambda^{2}+\epsilon_{\perp}^{2}}-f_{\mathrm{d}}$ is the
detuning,
$E_{\perp}^{\prime\prime}=\sqrt{(\alpha^{2}-\beta^{2})E_{x}^{2}+E_{y}^{2}}$ is
the electric fields perpendicular to the NV axis, and $E_{x}$ and $E_{y}$ are
the in-plane electric fields (See Supplementary materials for details). Here,
$\hat{L}_{\pm}$ are the rising (lowering) operators in
$\ket{\pm^{\prime}}_{\text{o}}$ basis. When $\Delta\sim 0$, the remained
second term in eq. (3) contributes to the Rabi splitting and the Rabi
oscillation, which will be discussed in the following.

Figure 3 (a) shows the experimental sequence. The microwave frequencies are
swept around frequency resonant to the transition between $\ket{0}$ and
$\ket{1}$, which is roughly estimated from the PLE measurement. At the same
time, the readout laser is applied to measure the population in $\ket{0}$. A
sufficiently long readout pulse increases the population in $\ket{1}$ through
the optical pumping as shown in the inset of Fig. 2 (b). When the microwave
frequency and the ground-state splitting are in resonance, the population in
$\ket{1}$ is transferred to $\ket{0}$, resulting in the increase of PL counts.
Figure 3 (b) shows the result of the ODER. The maximum of the spectrum is
12.84 GHz, which is in good agreement with the ground-state splitting obtained
from the PLE measurement. The full width at half maximum (FWHM) is 130 MHz,
which is due to the spectral diffusion in the ground and optically excited
states.

Using the resonance frequency (12.84 GHz) obtained from the ODER measurement,
we investigate the microwave power dependence of the Rabi splitting by
increasing the microwave power. Figure 3 (c) shows the color map of the PL
counts as functions of the relative frequency and the square root of the
microwave power. The maximum Rabi splitting is 2.6 GHz at 824 $\mu$W1/2. The
power dependence of the Rabi splitting directly corresponds to the power
dependence of the Rabi frequency, which is estimated to be 3.37 MHz/$\mu$W1/2
from the fitting. Also, the electric susceptibility are estimated to be
$\epsilon_{\perp}^{\mathrm{ac}}=$ 961 kHz/(V cm-1) (See Supplementary
materials for details). $\epsilon_{\perp}^{\mathrm{ac}}$ is larger than that
obtained using dc electric fields, 363 kHz/(V cm-1). The difference can be
attributed to the electric field screening in the dc measurement, which
decreased the effective electric fields at the NV center [17]. The value of
$\epsilon_{\perp}^{\mathrm{ac}}$ is in good agreement with that of the
optically excited state of NV-, 1.4 MHz/(V cm-1) [17].

Figure 4: The Rabi oscillation and the microwave power dependence of the Rabi
frequency. (a) Schematic of the experimental sequence for the Rabi
oscillation. The microwave pulse is applied between the first (initialization)
and second (readout) laser pulses. The microwave frequency is set to the
resonance frequency obtained from the ODER. The repetition times, $N$, are
100. (b) Normalized counts as a function of the microwave pulse width. The
orange dots are the data and the brown curve is the cosinusoidal fit. The Rabi
frequency is 87.8 MHz. (c) The Rabi frequency as a function of the square root
of the input power. The slope is 3.86 MHz/$\mu$W1/2.

### The Rabi oscillation and the Ramsey interference

Using the resonance frequency obtained from the ODER measurement, the Rabi
oscillations between $\ket{0}$ and $\ket{1}$ are observed. Figure 4 (a) shows
the experimental sequence for the measurement. The microwave pulse is applied
between the two laser pulses. Figure 4 (b) shows the Rabi oscillation at the
input microwave power of 828 $\mu$W. The Rabi frequency is 87.8 MHz. Compared
to the power required to drive the spin of NV- at a similar distance from the
electrode (several $\mu$m), the required power is three orders of magnitude
smaller. Figure 4 (c) shows the Rabi frequency as a function of the square
root of the microwave power. The Rabi frequency increases linearly with the
square root of the microwave power. The slope is 3.86 MHz/$\mu$W1/2, which is
in good agreement with the value obtained from the measurement of the Rabi
splitting, 3.37 MHz/$\mu$W1/2.

Figure 5: The Ramsey interference. (a) Schematic of the experimental sequence
for the Ramsey interference. The two microwave pulses are applied between the
first (initialization) and second (readout) laser pulses. The microwave width
for $\pi/2$ pulse is determined to be 4.6 ns from the measurement of the Rabi
oscillation. The repetition times, $N$, are 100. (b) The normalized counts as
a function of the delay time of the microwave pulse (free precession time).
The blue dots are the data and the blue curve is the fit with the function,
$A\mathrm{sin}(\omega t+\phi)\mathrm{exp}(-t/T_{2}^{*})+B$, where $A$, $B$,
and $\phi$ are some constants, $\omega$ is the angular frequency, and $T_{2}$
is the orbital coherence time. $T_{2}^{*}$ is 30.2 ns.

We also investigate the coherence time, $T_{2}^{*}$, of the orbital state
using the Ramsey interference. Figure 5 (a) shows the experimental sequence.
The two microwave $\pi/2$ pulses are applied between the initialize and
readout laser pulses. Figure 5 (b) shows the Ramsey interference with some
detuning as a function of free precession time. From the curve fit, the
detuning is 58 MHz and $T_{2}^{*}$ is 30.2 ns. Origins of the decoherence are
attributed to the surrounding nuclear spins [31, 32], the spectral diffusion
[33] in the ground state, and thermal phonons resonant with the orbital
transition as proposed for SiV [34]. The value of $T_{2}^{*}$ is similar to
the spin coherence time of SiV at around 5 K (30-50 ns) [35, 36, 18]. Thus, a
further increase of $T_{2}^{*}$ is expected by decreasing the temperature as
SiV and the dynamical decoupling.

## Discussion

In conclusion, we have measured the electric susceptibility of NV0 in the
ground state, which is comparable to that in the excited state of NV-. Also,
we have also succeeded in coherent control over the orbital state of NV0 using
electric fields. The required power for the Rabi oscillation is hundreds of
microwatts, which is three orders of magnitude smaller than that required for
the spin control using magnetic fields. The highly efficient control is
advantageous especially for the operation in a dilution refrigerator, opening
up the possibility for interfacing a superconducting qubit and a color center
via electric fields. If we can integrate a high impedance superconducting
microwave resonator with zero-point voltage fluctuation tens of microvolts
[37, 38], the single-photon coupling is expected to be tens of kilohertz,
enabling the microwave control using several hundred photons.

Moreover, the control using phonons is promising. Since the mode volume of the
phononic resonator can be designed to be several orders of magnitude smaller
than that of typical microwave resonators [29, 39], the single-phonon coupling
can exceed 1 MHz, enabling the realization of the strong coupling between a
single phonon and NV0.

Although the measured $T_{1}$ has been 137 ns at 5.5 K, it is expected to
extend to microseconds in the environment of a dilution refrigerator due to a
decrease in thermal phonons. If it reaches several microseconds, it is
sufficient to use for an interface between a superconducting quantum bit [30],
which is proposed for the silicon-vacancy center [29]. Despite a lot of
challenges posed by the spectral diffusion, integration of a microwave
resonator, and operation in a dilution refrigerator, our study demonstrates a
promising path toward hybrid quantum systems.

## Methods

### Sample fabrication

We use [100]-cut, electronic grade single-crystal diamond samples synthesized
by chemical vapor deposition (CVD) (Element Six). Before the fabrication of
the electrode, the diamond substrate is kept in a mixture of H2SO4 and HNO3 at
200∘C for 60 min to remove any surface contamination and terminate the surface
with oxygen. Then, Au (500 nm) / Ti (10nm) electrodes are formed on the
substrate by using photolithography processes.

### Low temperature confocal microscopy

All experiments are performed in a closed-cycle optical cryostat (Cryostation
s50, Montana Instruments) at 5.5 K under an ambient magnetic field. A copper
sample holder is mounted on XYZ piezoelectric nanopositioners (ANPx101$\times
2$, ANPz101) using a thermal link. The diamond sample is fixed onto a homemade
printed circuit board (PCB) using aluminum tapes. The electrodes on the
diamond are connected to the PCB via gold wire bonds. The PCB is attached to
the sample holder with screws.

Optical excitation and collection are done using a home-built confocal
microscope. An objective lens (LMPLFLN100X, Olympus), whose numerical aperture
is 0.8, is scanned over the sample using a XYZ piezoelectric nanopositioner
(P-517.3CD, Physik Instrumente).

A 515 nm green laser (Cobolt, Hubner Photonics) is used for off-resonant
excitation to search nitrogen-vacancy centers. A 575 nm yellow laser (DL-SHG
pro, Toptica Photonics) is used to resonantly excite NV0 for initialization
and readout. Typically, a 1 $\mu$s pulse is used and the power is 1 $\mu$W. A
637 nm red laser (DL pro, Toptica Photonics) is used to resonantly excite NV-
to convert the negatively charged state to the neutrally charged state. The
pulse width for the charge initialization is 100 $\mu$s and the power is 200
$\mu$W. Each laser joins to the main optical path using dichroic mirrors
(FF552-Di02-25x36: 515 nm, FF605-Di02-25x36: 575 nm, FF649-Di01-25x36: 637 nm,
Semrock) after purified through bandpass filters (FF02-520/28-25: 515 nm,
FF03-575/25-25: 575 nm, FF01-637/7-25: 637 nm, Semrock).

The power of 515 nm laser is controlled by setting the internal output power.
The powers of 575 nm and 637 nm lasers are controlled using variable fiber
optic attenuators (V450PA: 575 nm, V600PA: 637 nm, Thorlabs). The pulse
duration of 515 nm laser is directly set by the signal sent from a field
programmable gate array (FPGA). The pulse durations of 575 nm and 637 nm laser
are gated using acousto-optic modulators (SFO4903-S-M200 0.4C2C-3-F2P-01: 575
nm, SFO3916-S-M200-0.4C2E-3-F2P-02: 637 nm, Gooch and Housego) controlled by
the FPGA (PXIe-7820R, National Instruments).

The phonon sideband of NV0 and NV- is collected and measured using a single-
photon counter (SPCM-AQRH-14-FC, Excelitas). For the time-resolved
photoluminescence measurements, the signal from the single-photon counter is
sent to a time-resolved single-photon counter (Picoharp300, PicoQuant). The
time-bin is set to be 512 ps.

### DC and ac electronics for measurements

The dc voltages are generated using an arbitrary waveform generator (M3202A,
Keysight). The voltages are amplified using amplifiers (EVAL-ADHV4702-1CPZ,
Analog Devices). The microwave voltages up to 16 GHz ($>$ 25 GHz for a single
channel) are generated using another arbitrary waveform generator (M8195A,
Keysight). The microwave voltages are amplified using an amplifier
(ZVE-3W-183+, Mini Circuits). The microwave power at the sample is estimated
using a network analyzer (P9373A, Keysight). The loss inside the cryostat is
estimated from the reflection coefficient, $S_{11}$, inside the cryostat.
Assuming equal energy loss on the way to the sample and on the way back, half
of $S_{11}$ corresponds to the loss between the input port and the sample.

## Acknowledgements

This work was supported by a Japan Science and Technology Agency (JST)
Moonshot R$\&$D grant (JPMJMS2062) and by a JST CREST grant (JPMJCR1773). We
also acknowledge the Ministry of Internal Affairs and Communications (MIC) for
funding, research and development for construction of a global quantum
cryptography network (JPMI00316) and the Japan Society for the Promotion of
Science (JSPS) Grants-in-Aid for Scientific Research (20H05661, 20K20441).

## Authorship contribution

H. Kurokawa and Y. S. designed the experiments. H. Kurokawa wrote the
manuscript. K. W. and S. N performed the experiments. T. M and H. Kato
fabricated the electrical circuit. H. Kosaka supervised the project.

## References

  * Bernien  _et al._ [2013] H. Bernien, B. Hensen, W. Pfaff, G. Koolstra, M. S. Blok, L. Robledo, T. H. Taminiau, M. Markham, D. J. Twitchen, L. Childress, and R. Hanson, Nature 497, 86 (2013).
  * Pompili  _et al._ [2021] M. Pompili, S. L. N. Hermans, S. Baier, H. K. C. Beukers, P. C. Humphreys, R. N. Schouten, R. F. L. Vermeulen, M. J. Tiggelman, L. dos Santos Martins, B. Dirkse, S. Wehner, and R. Hanson, Science 372, 259 (2021).
  * Bhaskar  _et al._ [2020] M. K. Bhaskar, R. Riedinger, B. Machielse, D. S. Levonian, C. T. Nguyen, E. N. Knall, H. Park, D. Englund, M. Lončar, D. D. Sukachev, and M. D. Lukin, Nature 580, 60 (2020).
  * Nemoto  _et al._ [2014] K. Nemoto, M. Trupke, S. J. Devitt, A. M. Stephens, B. Scharfenberger, K. Buczak, T. Nöbauer, M. S. Everitt, J. Schmiedmayer, and W. J. Munro, Physical Review X 4, 031022 (2014).
  * Abobeih  _et al._ [2019] M. H. Abobeih, J. Randall, C. E. Bradley, H. P. Bartling, M. A. Bakker, M. J. Degen, M. Markham, D. J. Twitchen, and T. H. Taminiau, Nature 576, 411 (2019).
  * Dolde  _et al._ [2011] F. Dolde, H. Fedder, M. W. Doherty, T. Nöbauer, F. Rempp, G. Balasubramanian, T. Wolf, F. Reinhard, L. C. L. Hollenberg, F. Jelezko, and J. Wrachtrup, Nature Physics 7, 459 (2011).
  * Degen  _et al._ [2017] C. L. Degen, F. Reinhard, and P. Cappellaro, Reviews of Modern Physics 89, 035002 (2017).
  * 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).
  * Bar-Gill  _et al._ [2013] N. Bar-Gill, L. M. Pham, A. Jarmola, D. Budker, and R. L. Walsworth, Nature Communications 4, 1743 (2013).
  * Abobeih  _et al._ [2018] M. H. Abobeih, J. Cramer, M. A. Bakker, N. Kalb, M. Markham, D. J. Twitchen, and T. H. Taminiau, Nature Communications 9, 2552 (2018).
  * Abobeih  _et al._ [2022] M. H. Abobeih, Y. Wang, J. Randall, S. J. H. Loenen, C. E. Bradley, M. Markham, D. J. Twitchen, B. M. Terhal, and T. H. Taminiau, Nature 606, 884 (2022).
  * Nguyen  _et al._ [2019a] C. T. Nguyen, D. D. Sukachev, M. K. Bhaskar, B. Machielse, D. S. Levonian, E. N. Knall, P. Stroganov, R. Riedinger, H. Park, M. Lončar, and M. D. Lukin, Physical Review Letters 123, 183602 (2019a).
  * Tamarat  _et al._ [2008] P. Tamarat, N. B. Manson, J. P. Harrison, R. L. McMurtrie, A. Nizovtsev, C. Santori, R. G. Beausoleil, P. Neumann, T. Gaebel, F. Jelezko, P. Hemmer, and J. Wrachtrup, New Journal of Physics 10, 045004 (2008).
  * Maze  _et al._ [2011] J. R. Maze, A. Gali, E. Togan, Y. Chu, A. Trifonov, E. Kaxiras, and M. D. Lukin, New Journal of Physics 13, 025025 (2011).
  * Meesala  _et al._ [2018] S. Meesala, Y.-I. Sohn, B. Pingault, L. Shao, H. A. Atikian, J. Holzgrafe, M. Gündogan, C. Stavrakas, A. Sipahigil, C. Chia, R. Evans, M. J. Burek, M. Zhang, L. Wu, J. L. Pacheco, J. Abraham, E. Bielejec, M. D. Lukin, M. Atatüre, and M. Lončar, Physical Review B 97, 205444 (2018).
  * Chen  _et al._ [2018] H. Y. Chen, E. R. MacQuarrie, and G. D. Fuchs, Physical Review Letters 120, 167401 (2018).
  * Block  _et al._ [2021] M. Block, B. Kobrin, A. Jarmola, S. Hsieh, C. Zu, N. L. Figueroa, V. M. Acosta, J. Minguzzi, J. R. Maze, D. Budker, and N. Y. Yao, Physical Review Applied 16, 024024 (2021).
  * Maity  _et al._ [2020] S. Maity, L. Shao, S. Bogdanović, S. Meesala, Y.-I. Sohn, N. Sinclair, B. Pingault, M. Chalupnik, C. Chia, L. Zheng, K. Lai, and M. Lončar, Nature Communications 11, 193 (2020).
  * Batalov  _et al._ [2008] A. Batalov, C. Zierl, T. Gaebel, P. Neumann, I.-Y. Chan, G. Balasubramanian, P. R. Hemmer, F. Jelezko, and J. Wrachtrup, Physical Review Letters 100, 077401 (2008).
  * Wolfowicz  _et al._ [2020] G. Wolfowicz, F. J. Heremans, C. P. Anderson, S. Kanai, H. Seo, A. Gali, G. Galli, and D. D. Awschalom, Nature Reviews Materials 6, 906 (2020).
  * Manson and Harrison [2005] N. B. Manson and J. P. Harrison, Diamond and Related Materials 14, 1705 (2005).
  * Hauf  _et al._ [2011] M. V. Hauf, B. Grotz, B. Naydenov, M. Dankerl, S. Pezzagna, J. Meijer, F. Jelezko, J. Wrachtrup, M. Stutzmann, F. Reinhard, and J. A. Garrido, Physical Review B 83, 081304 (2011).
  * Grotz  _et al._ [2012] B. Grotz, M. V. Hauf, M. Dankerl, B. Naydenov, S. Pezzagna, J. Meijer, F. Jelezko, J. Wrachtrup, M. Stutzmann, F. Reinhard, and J. A. Garrido, Nature Communications 3, 729 (2012).
  * Doi  _et al._ [2014] Y. Doi, T. Makino, H. Kato, D. Takeuchi, M. Ogura, H. Okushi, H. Morishita, T. Tashima, S. Miwa, S. Yamasaki, P. Neumann, J. Wrachtrup, Y. Suzuki, and N. Mizuochi, Physical Review X 4, 011057 (2014).
  * Barson  _et al._ [2019] M. S. Barson, E. Krausz, N. B. Manson, and M. W. Doherty, Nanophotonics 8, 1985 (2019).
  * Baier  _et al._ [2020] S. Baier, C. E. Bradley, T. Middelburg, V. V. Dobrovitski, T. H. Taminiau, and R. Hanson, Physical Review Letters 125, 193601 (2020).
  * Zhang  _et al._ [2021] Q. Zhang, Y. Guo, W. Ji, M. Wang, J. Yin, F. Kong, Y. Lin, C. Yin, F. Shi, Y. Wang, and J. Du, Nature Communications 12, 1529 (2021).
  * Aslam  _et al._ [2013] N. Aslam, G. Waldherr, P. Neumann, F. Jelezko, and J. Wrachtrup, New Journal of Physics 15, 013064 (2013).
  * Neuman  _et al._ [2021] T. Neuman, M. Eichenfield, M. E. Trusheim, L. Hackett, P. Narang, and D. Englund, npj Quantum Information 7, 121 (2021).
  * Kurokawa  _et al._ [2022] H. Kurokawa, M. Yamamoto, Y. Sekiguchi, and H. Kosaka, Physical Review Applied 18, 064039 (2022).
  * Becker  _et al._ [2018] J. N. Becker, B. Pingault, D. Groß, M. Gündogan, N. Kukharchyk, M. Markham, A. Edmonds, M. Atatüre, P. Bushev, and C. Becher, Physical Review Letters 120, 053603 (2018).
  * Nguyen  _et al._ [2019b] C. T. Nguyen, D. D. Sukachev, M. K. Bhaskar, B. Machielse, D. S. Levonian, E. N. Knall, P. Stroganov, C. Chia, M. J. Burek, R. Riedinger, H. Park, M. Lončar, and M. D. Lukin, Physical Review B 100, 165428 (2019b).
  * Robledo  _et al._ [2010] L. Robledo, H. Bernien, I. van Weperen, and R. Hanson, Physical Review Letters 105, 177403 (2010).
  * Jahnke  _et al._ [2015] K. D. Jahnke, A. Sipahigil, J. M. Binder, M. W. Doherty, M. Metsch, L. J. Rogers, N. B. Manson, M. D. Lukin, and F. Jelezko, New Journal of Physics 17, 043011 (2015).
  * Rogers  _et al._ [2014] L. J. Rogers, K. D. Jahnke, M. H. Metsch, A. Sipahigil, J. M. Binder, T. Teraji, H. Sumiya, J. Isoya, M. D. Lukin, P. Hemmer, and F. Jelezko, Physical Review Letters 113, 263602 (2014).
  * Pingault  _et al._ [2014] B. Pingault, J. N. Becker, C. H. H. Schulte, C. Arend, C. Hepp, T. Godde, A. I. Tartakovskii, M. Markham, C. Becher, and M. Atatüre, Physical Review Letters 113, 263601 (2014).
  * Samkharadze  _et al._ [2016] N. Samkharadze, A. Bruno, P. Scarlino, G. Zheng, D. P. DiVincenzo, L. DiCarlo, and L. M. K. Vandersypen, Physical Review Applied 5, 044004 (2016).
  * Niepce  _et al._ [2019] D. Niepce, J. Burnett, and J. Bylander, Physical Review Applied 11, 044014 (2019).
  * [39] B. Kim, H. Kurokawa, H. Kosaka, and M. Nomura, arXiv:2305.08306 .