\section{Introduction}
The negatively charged nitrogen-vacancy (NV$^-$ or just NV) center in diamond has opened up a unique opportunity
for nanoscale nuclear magnetic resonance (NMR) spectroscopy and imaging.~\cite{MKS+13,SSP+13,MKC+14,HSR+15,DPL+15,LSU+16,APN+17}
NV's performance as a room temperature quantum sensor for nanoNMR is crucially dependent on
its depth from the diamond surface ($d_{\mathrm{NV}}$), as the NV center closer to the surface exhibits a stronger interaction with NMR analyte placed on the surface.
Furthermore, near-surface NV centers lead to efficient couplings with light,
which also adds impetus to other applications such as nanophotonics, plasmonics, and quantum communication.~\cite{TMO+13,SMZ+16,AEV+18,AHWZ18,WEH18,BGF+19}

There are several approaches to creating near-surface NV centers.
Most commonly used is low-energy N$^+$ ion implantation, typically with the acceleration energy less than 5~keV.~\cite{PNJ+10,OPC+12,dOAW+17,FBH+18,SDS+19}
Other methods such as nitrogen $\delta$-doping by chemical vapor deposition (CVD) and precise surface etching by plasma or thermal oxidation have also been demonstrated.~\cite{OHB+12,ORW+13,MDD+14,LPMD14,KMS+14,CGO+15,dOMW+15,ZZW+17}
Of particular interest to these platforms is the relation between $d_{\mathrm{NV}}$ and $T_2$, the coherence time of NV's electronic spin.
Longer $T_2$ is required to achieve better magnetic sensitivity, whereas near-surface NV centers tend to suffer short $T_2$.

In Ref.~\onlinecite{ISS+17}, some of the present authors reported that
near-surface NV centers can also be created by comparatively high energy (10~keV) N$^+$ ion implantation,
when combined with a SiO$_2$ screening mask deposited on the diamond surface.
The defining feature of this method is that the region of the highest N$^+$ density is located at the surface.
We showed this by Monte Carlo simulations, and experimentally confirmed the creation of NV centers
(see Sec.~\ref{sec_sample} for a more detailed description of our method and the result of Monte Carlo simulations).
The resulting N$^+$ distribution in our method is in contrast with the case for the standard (maskless) low-energy ion implantation,
in which the implanted N$^+$ ions are approximately Gaussian-distributed inside of the diamond.
An experimental confirmation of the distribution of the NV centers along the depth direction (depth profile),
as well as the relation between $T_2$ and $d_{\mathrm{NV}}$, are highly desired but have not yet been reported so far.
The present work addresses these issues.
Specifically, we focus on single NV centers and determine $d_{\mathrm{NV}}$ from NMR of proton ensembles on the surface, with direct application to nanoNMR in mind.~\cite{AS18,SIA18}

The rest of this paper is organized as follows.
In Sec.~\ref{sec_density}, we determine the density and yield of NV centers created by our method.
Section~\ref{sec_t2} presents the main result, the relation between $T_2$ and $d_{\mathrm{NV}}$,
and we analyze the depth profile of the NV centers by assembling the data from as many as 141 NV centers.
In Sec.~\ref{sec_noise}, we perform noise spectroscopy to investigate the source of the noise that limits $T_2$ in our sample.
The diamond surface is examined by atomic force microscopy (AFM) and X-ray photoelectron spectroscopy (XPS) in Sec.~\ref{sec_surface}.
We conclude in Sec.~\ref{sec_conclusion}.
To avoid digressing from the main subject, technical aspects and supporting data are presented in Appendices for interested readers.
They include the sample preparation procedure and the Monte Carlo simulations of N$^+$ ion implantation with screening masks (Sec.~\ref{sec_sample}), 
our experimental setup and microwave pulse sequences (Sec.~\ref{sec_pulse}),
and the result from double electron--electron resonance (DEER) measurements (Sec.~\ref{sec_deer}).
Although many of them are now becoming common tools among the researchers working on NV-based quantum sensing, it will still be beneficial to provide concise summaries.
We also refer to Ref.~\onlinecite{MSI+20} as a resource for further details about our experimental setup and techniques extensively used in the present work.

\section{Density and yield\label{sec_density}}
Figures~\ref{fig1}(a--d) show fluorescence images (20~$\times$~20~$\mu$m$^2$) of diamond surfaces,
where the SiO$_2$ layers of thickness $t$ = 52.3, 57.6, 64.1, and 69.1~nm were deposited prior to N$^+$ ion implantation.
\begin{figure}
\begin{center}
\includegraphics{fig1.eps}
\caption{(a--d) Fluorescence images (20~$\times$~20~$\mu$m$^2$) of diamond surfaces under the laser power of 0.6~mW.
$t$: thickness of the SiO$_2$ layer.
kcps: kilo-counts per second.
(e) NV density as a function of $t$.
(f) Yield as a function of $t$.
\label{fig1}}
\end{center}
\end{figure}
It should be noted that these SiO$_2$ layers had been removed when these images were taken (see Sec.~\ref{sec_sample}).
Each image shows multiple bright spots, the number of which is decreasing as $t$ increases.
To confirm the presence of single NV centers, continuous wave optically detected magnetic resonance (CW ODMR) was performed on individual bright spots.
For the areas with $t$ = 64.1 and 69.1~nm, all the bright spots in the searched area were examined and the number of ODMR-active points were counted.
The areal density was then calculated by dividing it by the searched area.
For the $t$ = 52.3 and 57.6~nm areas, where many more bright spots exist, we set a threshold count rate, above which we judged as an emission from a single NV center,
and the number of above-threshold spots were counted by a computer program.
The threshold count rates were 2.0~$\times$~10$^4$~kcps for $t$ = 52.3~nm and 1.7~$\times$10$^4$~kcps for $t$ = 57.6~nm,
which were the minimum photon count rates among a subset of the bright spots on which CW ODMR was performed and the signals were observed.
The result is summarized in Fig.~\ref{fig1}(e).

We also calculate the yield (efficiency of conversion from N$^+$ ions into NV$^-$ centers) as shown in Fig.~\ref{fig1}(f).
It is given by the ratio of the areal density divided by the implanted N$^+$ density deduced from Monte Carlo simulations.
The low yield of 0.1--0.4\% is consistent with our previous observation,~\cite{ISS+17}
and similar to the values obtained in low energy ion implantation ($\approx$ 3~keV).~\cite{PNJ+10}
The cause of the latter case is attributed to the lack of vacancies to pair up with nitrogen atoms to form NV centers,
but it has been shown that the yield can be improved by carefully adjusting annealing, chemical or plasma etching conditions.~\cite{OSP+13,YUW+13,AHA+14}
This aspect has not been optimized yet in our case, and can be improved in the future.

\section{Coherence time and depth\label{sec_t2}}
Figure~\ref{fig2}(a) shows an example of Hahn echo decay [see Fig.~\ref{fig7}(a) for the pulse sequence].~\cite{H50}
\begin{figure}
\begin{center}
\includegraphics{fig2.eps}
\caption{(a) Hahn echo decay.
(b) NMR spectrum.
(c) Relation between $T_{2,\mathrm{echo}}$ and $d_{\mathrm{NV}}$.
The dashed line is the measurement limit as determined by Eq.~(\ref{eq_t2limit}), whereas the dash--dot line indicates $T_2$ = 3~$\mu$s.
They define the region (unpainted) that is ``accessible'' in our experiments.
See the main text for detail.
The filled circle indicates the NV center with which the Hahn echo decay in (a), the NMR spectrum in (b), and the noise spectroscopy data in Fig.~\ref{fig4} were taken.
The filled triangle indicates the NV center with which the DEER data in Fig.~\ref{fig8} were taken.
\label{fig2}}
\end{center}
\end{figure}
The static magnetic field $B_0$ = 23.2~mT was applied parallel to the NV symmetry axis,
and the microwave frequency was set at $\omega_{0,-1}/2\pi$ = 2218.2~MHz, corresponding to the $m_S = 0 \leftrightarrow -1$ transition.
$P_0$ is the probability of the final spin state being $m_S$ = 0, which we calculate from calibrated photon counts.~\cite{MSI+20}
We describe a general coherence decay curve as
\begin{equation}
P_0 = \frac{1}{2} + \frac{ \mathcal{C} }{2} \exp \left[ -\left( \frac{ t_{\mathrm{tot}} }{T_2} \right)^p \right].
\label{eq_t2}
\end{equation}
The form of $\mathcal{C}$ depends on a type of measurements, $t_{\mathrm{tot}}$ is the total time of the pulse sequence applied.
Fitting parameters are $T_2$ and $p$ (stretched exponent).
In the case of Hahn echo, $\mathcal{C}_{\mathrm{echo}}$ = 1 and $t_{\mathrm{tot}}$ = 2$\tau$.
The solid line in Fig.~\ref{fig2}(a) is the fit, giving $T_{2,\mathrm{echo}}$ = 27.1~$\mu$s and $p$ = 0.98.

We then apply a multipulse sequence to determine $d_{\mathrm{NV}}$ from NMR of protons on the surface.
We use XY4, XY8, and XY16 sequences (see Sec.~\ref{sec_pulse} for the definitions).~\cite{GBC90}
When the total number of $\pi$ pulses in a sequence is $N$, it is denoted as XY$k$-$N$ ($k$ = 4, 8, or 16 and $N$ is an integer multiple of $k$).
The decay is given by setting $t_{\mathrm{tot}}$ = $N\tau$, and the NMR spectrum is described by~\cite{PDC+16} 
\begin{eqnarray}
& & \mathcal{C}_{\mathrm{NMR}} = \\
& & \exp \left\{ - 2 \left( \frac{ \gamma_{\mathrm{e}} B_{\mathrm{rms}} N\tau }{ \pi } \right)^2 
\mathrm{sinc}^2 \left[ \frac{N \tau}{2} \left( \omega_{\mathrm{n}} - \frac{\pi}{\tau} \right) \right] \right\}. \nonumber
\end{eqnarray}
Here, $\gamma_{\mathrm{e}}/2\pi$ = 28~MHz/mT is the gyromagnetic ratio of the NV spin, and
$\omega_{\mathrm{n}}/2\pi = \gamma_{\mathrm{h}} B_0/2\pi$ is the proton NMR frequency
with $\gamma_{\mathrm{h}}/2\pi$ = 42.577~kHz/mT the gyromagnetic ratio of the $^1$H nuclear spin.
$B_{\mathrm{rms}}$ is the root-mean-square nuclear dipolar magnetic field integrated over the semi-infinite volume (i.e., outside of the diamond),
and can be calculated analytically as~\cite{SSP+13,LPMD14,RDO+14,PDC+16} 
\begin{equation}
B_{\mathrm{rms}} = \frac{\mu_0 \hbar \gamma_{\mathrm{h}} }{4 \pi} \sqrt{ \frac{ 5 \pi \rho}{ 96 \, d_{\mathrm{NV}}^3 } }
\end{equation}
with $\mu_0$ the permeability of vacuum, $\hbar$ the reduced Planck constant, and
$\rho$ the proton density of the immersion oil on the surface, known to be 6~$\times$~10$^{28}$~m$^{-3}$ in the present case (Olympus IMMOIL-F30CC).

Figure~\ref{fig2}(b) is the NMR spectrum, taken with the NV center of Fig.~\ref{fig2}(a) and the XY16-128 sequence.
The stretched exponential decay component is subtracted (normalized $P_0$, or $P_{0,\mathrm{norm}}$),
and the horizontal axis is $(2\tau)^{-1}$ in the frequency unit.
The NMR dip appears at 0.988~MHz, agreeing with the calculated $\omega_{\mathrm{n}}/2\pi$. 
From the fit, $d_{\mathrm{NV}}$ is determined to be 17.4~nm.
(Consult Fig.~10 of Ref.~\onlinecite{MSI+20} for another example of the NMR spectrum measured with a different, shallower NV center in the same sample.)

We repeat the measurements on different NV centers, and the relation between $T_{2,\mathrm{echo}}$ and $d_{\mathrm{NV}}$,
obtained from 43 NV centers in total, is plotted in Fig.~\ref{fig2}(c).
[The point with the filled circle corresponds to the NV center examined in Figs.~\ref{fig2}(a,b).]
We note again that, although the NV centers are grouped by $t$ (SiO$_2$ thickness), the SiO$_2$ layer had been removed when these measurements were performed.
A general trend is clearly discerned: the deeper the NV center, the longer the coherence,
as observed in previous works using different methods to create NV centers.~\cite{dOAW+17,FBH+18,SDS+19}
For instance, we find the observed $T_{2,\mathrm{echo}}$--$d_{\mathrm{NV}}$ relation is quantitatively close to that obtained in Ref.~\onlinecite{SDS+19} (without high temperature annealing).
It must be emphasized, however, that in our case a large fraction of NV centers exhibited short $T_2$ insufficient to observe clear NMR spectra.
Consequently, we could not determine their depths, and such NV centers are not shown in Fig.~\ref{fig2}(c).
We will include them when estimating the depth profile,
but before doing so we discuss an experimental condition to detect NMR and set the region in the $T_{2,\mathrm{echo}}$--$d_{\mathrm{NV}}$ space
where the experimental observation can be made.

The condition is simply that the signal $\mathcal{S}$ be larger than the noise $\mathcal{N}$.
The NMR signal arises from the component of $P_0$ after subtracting the coherence decay,
therefore $\mathcal{S} = \frac{1}{2} (1-\mathcal{C}_{\mathrm{NMR}}) e^{-(N\tau/T_2)^p}$.
We set $p$ = 1 for simplicity (also consistent with experimental observation) and $\tau = \pi/\omega_{\mathrm{n}}$ = 506.2~ns, corresponding to the resonance condition.
Experimentally, the number of $\pi$ pulses $N$ is determined from the balance between the signal strength and the coherence time,
and we usually found $N\tau/T_{2,\mathrm{echo}} \approx 2$ to be an appropriate operating condition.
For instance, the condition for Fig.~\ref{fig2}(b) is $N\tau/T_{2,\mathrm{echo}}$ = 128~$\times$~506.2~ns/27.1~$\mu$s $\approx$ 2.4.
Here, $T_{2,\mathrm{echo}}$ is used, as we do not know $T_2$ under dynamical decoupling in advance.

The measurement noise $\mathcal{N}$ arises from the Poisson distribution of the photon counts $n$, the variance of which is given by $\sqrt{n}$.
The signal photon count $n_{\mathrm{s}}$ is defined as $n_{\mathrm{s}} = n_1 + (n_0 - n_1) P_0$,
where $n_0$ and $n_1$ are the photon counts from the $m_S$ = 0 and $-$1 states, respectively, and are mutually related by $n_1 = (1-c) n_0$, with $c$ the contrast.
We obtain $P_0 = (n_{\mathrm{s}} - n_1)/(c n_0) $
and the fluctuation of $n_{\mathrm{s}}$, given by $\sqrt{n_{\mathrm{s}}}$, leads to $\mathcal{N} = \sqrt{ n_{\mathrm{s}} }/(c n_0) \approx (c \sqrt{n_0})^{-1}$.
Combining these, we write the condition $\mathcal{S}/\mathcal{N} = 1$ as
\begin{equation}
T_{2,\mathrm{echo}} = \frac{ \pi }{ 2\sqrt{2} \gamma_e B_{\mathrm{rms}} } \sqrt{-\ln \left( 1- \frac{2 e^2}{ c \sqrt{n_0} } \right) }.
\label{eq_t2limit}
\end{equation}
Typical values for $c$ and $n_0$ in the present experimental condition are $c \approx$ 0.2 and $n_0 \approx 5 \times 10^5$ counts.
Substituting these values, we draw the dashed line in Fig.~\ref{fig2}(c), below which the depth cannot be determined.
In addition to this constraint, short $T_{2,\mathrm{echo}}$ limits the number of pulses we can apply.
For instance, to satisfy $N\tau/T_{2,\mathrm{echo}} \approx 2$ with $N$ = 8, $T_{2,\mathrm{echo}}$ needs to be longer than 2~$\mu$s.
In fact, from NV centers with $T_{2,\mathrm{echo}}$ shorter than 3~$\mu$s we were not able to observe clear NMR signals
(and the extension of $T_2$ by dynamical decoupling is often not effective),
even though they are likely to be located close to the surface.
This motivates us to draw another bound $T_{2,\mathrm{echo}}$ = 3~$\mu$s in Fig.~\ref{fig2}(c).
Although the argument here is manifestly crude, the defined ``accessible'' region (unpainted) reasonably coincides with the experimental observation.

From the dataset in Fig.~\ref{fig2}(c), it is straightforward to produce histograms
showing the number of NV centers (occurrence) as a function of $d_{\mathrm{NV}}$ [Figs.~\ref{fig3}(a--d)].
\begin{figure}
\begin{center}
\includegraphics{fig3.eps}
\caption{(a--d) Histograms showing the distributions of the 43 {\it depth-determined} NV centers.
(e--h) NV densities estimated from the experiments ($\bigcirc$), 98 {\it depth-undetermined} NV centers included, and the simulations (solid lines).
See the main text for detail.
\label{fig3}}
\end{center}
\end{figure}
Ideally, these histograms are to be directly compared with the depth profile from Monte Carlo simulations (Fig.~\ref{fig6}).
However, as mentioned above, the fact that the data points in Fig.~\ref{fig2}(c) are limited to a certain region suggests
that they alone are not sufficient to deduce a firm conclusion about the depth profile.
There are indeed another 98 NV centers we measured but $d_{\mathrm{NV}}$ were not determined.
This number breaks down into 44 ($t$ = 52.3~nm), 33 (57.6~nm), 10 (64.1~nm), and 11 (69.1~nm),
not negligible compared with the numbers of the depth-determined NV centers, 13 (52.3~nm), 12 (57.6~nm), 9 (64.1~nm), and 9 (69.1~nm).

To include the former group in the histograms, we assume, judging from the trend observed in Fig.~\ref{fig2}(c), that they are located at $d_{\mathrm{NV}} \leq$ 5~nm.
In Figs.~\ref{fig3}(e--h), the solid lines are the Monte Carlo simulations in Fig.~\ref{fig6} multiplied by the respective yields in Fig.~\ref{fig1}(f).
The circle points ($\bigcirc$) for $d_{\mathrm{NV}} >$ 5~nm are the histograms multiplied by the respective NV densities.
Similarly, those for $d_{\mathrm{NV}} \leq$ 5~nm are based on the numbers of the NV centers whose depths left undetermined, but divided by five to reflect the fact that
we are agnostic on which ``bin'' (the width of 1~nm) they are in.
Thus, the horizontal error bar for them is set as $d_{\mathrm{NV}}$ = 2.5~$\pm$ 2.5~nm (the other error bars are set as the bin width).
From Figs.~\ref{fig3}(e--h), it is interpreted, though not conclusive, that the profile is concentrated toward the surface,
unlike the standard ion implantation without a screening mask.

\section{Noise spectroscopy\label{sec_noise}}
In our sample, it is observed that the NV centers with short $T_{2,\mathrm{echo}}$ also tend to be short-lived, degrading $T_{2,\mathrm{echo}}$ during the experiments.
In one case, initially $T_{2,\mathrm{echo}}$ was 6~$\mu$s, but it reduced continuously to 2~$\mu$s in 5 hours, then fluctuates between 1 and 2~$\mu$s.
Another case shows that initial $T_{2,\mathrm{echo}}$ of 18~$\mu$s slowly degrades to 12~$\mu$s in the time scale of 80~hours.
These are likely to arise from the coupling with the surface, and we next examine the origin of the noise that affects the NV spin.
For this purpose, we would ideally hope to use shallower NV centers, but as mentioned above, they suffer degradation of $T_{2,\mathrm{echo}}$,
making a clean set of experiments difficult.
Nonetheless, the trend in Fig.~\ref{fig2}(c) indicates that the plotted NV centers share the same noise source.
We therefore chose to use the NV center with longest $T_{2,\mathrm{echo}}$ found in the $t$ = 52.3~nm area [the filled circle in Fig.~\ref{fig2}(c)].
For the noise spectroscopy, $B_0$ was fixed at 23.2~mT, where $T_2$ and $d_{\mathrm{NV}}$ had also been measured.

We first perform dynamical decoupling by increasing $N$ up to 512.
The decay curves are shown in Fig.~\ref{fig4}(a).
\begin{figure}
\begin{center}
\includegraphics{fig4.eps}
\caption{(a) Coherence decay curves under dynamical decoupling.
The $N$ = 1 data were taken by the Hahn echo (the time axis is given by 2$\tau$),
the $N$ = 8 data by XY8, and the $N \geq$ 16 data by XY16-$N$ (the time axis is given by $N\tau$).
The inset: $T_2$ as a function of $N$.
The horizontal bars indicate the 95\% confidence limits.
(b) Relaxation signals for the SQ ($\triangle$) and DQ ($\square$) channels.
(c) Noise spectrum $S_{\mathrm{e}}(\omega)$ deduced from the spectral decomposition of the dynamical decoupling data ($\bigcirc$),
the SQ relaxation rate $\Omega$ ($\triangle$), and the DQ relaxation rate $\gamma$ ($\square$).
All the measurements were performed at $B_0$ = 23.2~mT.
\label{fig4}}
\end{center}
\end{figure}
Although dynamical decoupling improves the coherence to some extent, the effect is not significant.
$T_2^{(N)}$, $T_2$ by $N$-pulse dynamical decoupling, shows a scaling and saturating behavior described by
\begin{equation}
T_2^{(N)} = T_2^{(1)} [ N_{\mathrm{sat}}^s + (N^s - N_{\mathrm{sat}}^s) e^{-N/N_{\mathrm{sat}}} ],
\end{equation}
and the fit gives $N_{\mathrm{sat}}$ = 105.2, $s$ = 0.10, and $T_2^{(\infty)}$ = $T_2^{(1)} N_{\mathrm{sat}}^{s}$ = 43.0~$\mu$s [the inset of Fig.~\ref{fig4}(a)].
$T_2^{(1)}$ corresponds to $T_{2,\mathrm{echo}}$, but is treated as a parameter here.
The scaling exponent $s$ of 2/3 is expected for an environment that fluctuates slowly compared with $T_2$.~\cite{BPB+12,RMU+15}
The small $s \ll$ 2/3 as well as the saturated $T_2$ observed here indicate the presence of broadband noise.

This observation is reinforced by spectral decomposition.
Under a few simplifying assumptions, the spectral component at $\omega = \pi/\tau$ is given by~\cite{BPB+12,RMU+15}
\begin{equation}
S_{\mathrm{DD}}(\omega) = -\frac{ \pi \ln(2P_0 -1)}{ N \tau }.
\end{equation}
One can compute the noise spectrum directly from the dynamical decoupling data.
Figure~\ref{fig4}(c) shows the spectral function $S_{\mathrm{e}}(\omega)$,
related to $S_{\mathrm{DD}}(\omega)$ by $S_{\mathrm{e}}(\omega) = 2S_{\mathrm{DD}}(\omega)/d_{\parallel}^2$,
where $d_{\parallel}$ = 0.35~Hz~cm/V is the electric dipole moment parallel to the NV symmetry axis.~\cite{MAJ17}
The reason for this conversion is to enable a comparison between the (effective) electric and magnetic field noises.
Note that in Fig~\ref{fig4}(c) the data points close to $(2P_0 -1)$ = 0 or 1 that are insensitive to the change in $\tau$ have been removed.
Further analysis usually invokes a Lorentzian fit~\cite{BPB+12,RMU+15}
\begin{equation}
S_{\mathrm{DD}}(\omega) = \frac{ \langle E_{\perp}^2 \rangle \tau_{\mathrm{c}} }{ \pi [1+(\omega \tau_{\mathrm{c}})^2] } 
\end{equation}
or a sum of multiple Lorentzians, 
where $\langle E_{\perp}^2 \rangle$ is the average coupling strength of the environment to the NV center,
and $\tau_{\mathrm{c}}$ is the correlation time of the environment.
In the present case, the spectrum is almost flat (typically 1--2$\times$10$^{10}$~V$^2$/m$^2$Hz)
up to $\omega/2\pi \approx$ 10$^7$~Hz in the double logarithmic plot,
suggesting that the noise is white in the probed frequency range ($\omega \tau_{ \mathrm{c}} \ll 1$).
This makes fitting to a Lorentzian curve arbitrary, and we can only infer $\tau_{ \mathrm{c}}$ to be shorter than 100~ns.
Still, the short $\tau_{ \mathrm{c}} \ll T_2$ is consistent with the small $s$.

We also probe the noise by relaxation measurements through single quantum (SQ) and double quantum (DQ) channels
as shown in Fig.~\ref{fig4}(b) [see Fig.~\ref{fig7}(c) for the pulse sequence].~\cite{MAJ17}
The relaxation times of the two channels are almost the same,
$T_{1,\mathrm{SQ}}$ = 1.06~ms and $T_{1,\mathrm{DQ}}$ = 1.16~ms,
consistent with the previous observations at similar magnetic fields.~\cite{MAJ17,SDS+19,GCK20}
We can then deduce the SQ and DQ relaxation rates, $\Omega$ and $\gamma$,
using the relations $T_{1,\mathrm{SQ}} = (3 \Omega)^{-1}$ and $T_{1,\mathrm{DQ}} = (\Omega+2\gamma)^{-1}$.~\cite{MAJ17}
We obtain $\Omega$ = 316~Hz, corresponding to the magnetic noise probed at $\omega_{0,-1}/2\pi$ = 2218.2~MHz,
and $\gamma$ = 272~Hz, corresponding to the electric noise probed at $\omega_{-1,1}/2\pi$ = 1301.5~MHz (the separation of the $m_S$ = $-$1 and 1 states at 23.2~mT).
$\Omega$ and $\gamma$ are converted into the spectral components as
$2 \Omega/d_{\parallel}^2$ = 5.16$\times$10$^7$~V$^2$/m$^2$Hz and
$\gamma/d_{\perp}^2$ = 9.41$\times$10$^{3}$~V$^2$/m$^2$Hz, respectively,
with $d_{\perp}$ = 17~Hz~cm/V the electric dipole moment perpendicular to the NV symmetry axis.~\cite{MAJ17,SDS+19}
The combined plot in Fig.~\ref{fig4}(c) certifies that in the wide frequency region the dominant contribution comes from the magnetic noise,
which is presumably due to surface defects, as has been argued in the case of low energy ion implantations.~\cite{RMU+15,MAJ17,SDS+19}

\section{Surface condition\label{sec_surface}}
We have found that our method of creating near-surface NV centers provides the $T_2$--$d_{\mathrm{NV}}$ relation comparable to other existing methods.
Yet, the SiO$_2$ deposition, in which an amorphous SiO$_2$ source (melting point $\approx$ 1700$^{\circ}$C) is heated by an electron beam
and sublimated SiO$_2$ is directed to the diamond substrate in vacuum,
is the process unique to ours, and it is of interest to examine an effect it might have on the diamond surface.
Here, we perform AFM and XPS on two areas, one in which SiO$_2$ was not deposited and ion implantation was blocked (denoted as $t$ = 0),
and the other in which ion implantation was performed through SiO$_2$ layers (denoted as $t \neq$ 0).
These measurements were performed after all the ODMR-based experiments were done (including DEER presented in Sec.~\ref{sec_deer}),
and therefore the $t$ = 0 area was subject to the same chemical and annealing processes (described in Sec.~\ref{sec_sample}) as the $t \neq$ 0 area.

Figure~\ref{fig5}(a) compares AFM topographic images (500~$\times$~500~nm$^2$) taken at the $t$ = 0 and $t \neq$ 0 areas.
\begin{figure}
\begin{center}
\includegraphics{fig5.eps}
\caption{(a) AFM topographic images (500~$\times$~500~nm$^2$) of the diamond surfaces.
(b) The normalized XPS spectra showing the carbon 1$s$ (left) and oxygen 1$s$ (right) peaks.
The spectra for $t \neq$ 0 are shifted by 0.1 for clarity. 
\label{fig5}}
\end{center}
\end{figure}
To be more precise, the latter belongs to the $t$ = 52.3~nm area, but is not correlated with the area on which the ODMR experiments were carried out.
In these particular areas shown in Fig.~\ref{fig5}(a), the root-mean-square roughness, $R_{\mathrm{rms}}$, is 0.11~nm for $t$ = 0 and 0.14~nm for $t \neq$ 0.
The measurements were repeated on different locations, and $R_{\mathrm{rms}}$ averaged over the measured areas were 0.17~nm for both $t$ = 0 and $t \neq$ 0,
indicating that the SiO$_2$ deposition and subsequent ion implantation have not altered the surface condition significantly.
Nonetheless, the change in the surface morphology, with more `wrinkles' in the latter, is discernible.
It is not clear yet whether these wrinkles affect the coherence properties of the NV spins.

Figure~\ref{fig5}(b) shows XPS spectra of carbon 1$s$ and oxygen 1$s$.
The C 1$s$ peak at 285~eV arises primarily from $sp^3$ (in bulk),
and the asymmetric shape indicates the presence of other bonds such as
C--O, C=O ($>$ 285~eV), and $sp^2$ ($<$ 285~eV).~\cite{BDM+14,SDS+19}
The O 1$s$ peak at 532~eV has a broad linewidth, presumably composed of 
C--O--C, C--C--H, and C=O ($<$ 532~eV).~\cite{BDM+14,SDS+19}
Despite the limited spectral resolution, we can conclude that the chemical composition of the surface is essentially the same between the $t$ = 0 and $t \neq$ 0 areas.

\section{Conclusion and outlook\label{sec_conclusion}}
To conclude, we created single nitrogen-vacancy (NV) centers by 10-keV N$^+$ ion implantation into diamond via SiO$_2$ layers of $t$ = 52.3, 57.6, 64.1, and 69.1~nm,
and examined the relation between $T_{2,\mathrm{echo}}$ and $d_{\mathrm{NV}}$.
The depth profiles incorporating the NV centers with $T_2 <$ 3~$\mu$s were found to be consistent with Monte Carlo simulations.
The effect of the surface on $T_2$ was examined by a combination of noise spectroscopy, surface topography, and XPS.

Although our method is found to give the $T_2$--$d_{\mathrm{NV}}$ relation comparable to other existing methods,
we also conceive that the surface undergoes a slight change in morphology due to SiO$_2$ deposition and $T_2$ becomes unstable and/or is shortened over time.
Recently, these aspects have been investigated in great detail by Sangtawesin {\it et al.},~\cite{SDS+19}
who emphasized the importance of keeping the surface purity and the quality of surface morphology in each step of sample preparation.
Longer coherence was achieved by combining annealing at 800$^{\circ}$C, at 1200$^{\circ}$C, and in oxygen atmosphere,
which are thought to remove residual vacancy complexes due to ion implantation and lead to highly ordered oxygen terminated surface.
While we did not adopt annealing at 1200$^{\circ}$C, this additional process is reported to improve $T_2$ several times longer (e.g., a few $\mu$s to ten $\mu$s)
and achieve $R_{\mathrm{rms}}$ of 0.06~nm (compared with our $R_{\mathrm{rms}}$ = 0.17~nm).
Given that without 1200$^{\circ}$C the $T_2$--$d_{\mathrm{NV}}$ relation in the present work and that in Ref.~\onlinecite{SDS+19} are similar,
we expect the process to provide an important clue for stabilizing and extending $T_2$ of NV center within 5~nm from the surface.

\section*{Acknowledgement}
K.S. was supported by the JSPS Grant-in-Aid for Research Fellowship for Young Scientists (DC1), Grant No.~JP17J05890.
T.T. was supported by the JSPS Grant-in-Aid for Scientific Research (KAKENHI) (B) Grant No.~20H02187 and (B) Grant No.~19H02617,
(S) Grant No.~16H06326
the JST CREST (JPMJCR1773), and MEXT Q‐LEAP (JPMXS0118068379).
K.M.I. was supported by the JSPS KAKENHI (S) Grant No.~26220602 and (B) Grant No.~19H02547,
the JST Development of Systems and Technologies for Advanced Measurement and Analysis (SENTAN),
and the Spintronics Research Network of Japan (Spin-RNJ).
The authors thank the technical staff of Central Service Facilities for Research at Keio University
for their assistance in XPS measurements.
The data that support the findings of the present work are available from the corresponding author upon reasonable request.

