NV centres by vacancies trapping in irradiated diamond: experiments and modelling

The negatively charged nitrogen-vacancy $\mathrm{\left(NV^{-} \right)}$ centre in diamond is a nanoscale quantum sensor [1] with a wide range of applications comprising three-dimensional current density imaging [2], nuclear magnetic resonance [3], molecular imaging [4] and condensed matter physics [5]. Compared to other quantum sensors, a device based on NV centre ensembles in diamond has the advantage of operating over a wide range of temperatures with a record-breaking sensitivity to volume ratio [6]. The sensitivity of a NV sensor strongly depends on spin density used for quantum sensing and on their proximity to diamond surface. Maximising this density without compromising decoherence time $\mathrm{(T_{2})}$ or signal contrast is a key challenge in the field [7].

One promising method, based on in-situ doping during CVD growth on (111)-oriented substrates, allows to achieve a preferred alignment of NVs with accuracy ${\sim}94\%-97\%$ [8, 9]. Although the NV centres produced by this method have shown an enhanced sensitivity, its major limit is related to a low conversion efficiency without additional radiation exposure for NV activation [10].

Nitrogen implantation followed by an annealing process is one of the most commonly used methods to generate a higher density of NVs in diamond with a nanometric precision [11–15]. Besides its spatial resolution, low energy implantation (few keV), required for near-surface NV positioning, shows a low NV yield. The yield of this method is in a range of few percent for low energy implantation, which is relevant for sensing applications. A higher yield to almost $50\%$ can be achieved by co-implantation of donor impurities, which mediate the charge state of defects and their trapping efficiency at thermal annealing [16]. Still, ion implantation technique suffers from residual lattice damage along the ion track. This can create di-vacancies or/and vacancy chains during a following annealing step, thus leading to a reduction of the spin coherence time [17, 18].

The properties of near-surface NVs by ion implantation can be improved to some extend using isotopically purified materials [11, 15, 19–21] and surface manipulation based on plasma [22, 23] and acid boiling [15, 24] reducing the effect of surface damage for NVs spin degradation.

Recently, a combination of CVD growth of N-doped layers and diffusion of vacancies from a remote source by ion irradiation (He, C) has been considered a valid alternative to nitrogen implantation [25–27]. A combination of nitrogen and helium implantation has been also applied for improving the properties of shallow NV ensembles [28]. In particular, the ion irradiation technique represents an attempt to separate the location of radiation damage from the NVs induced by diffusion of vacancies. In most of the experiments, however, the P1-to-NV conversion efficiency by such methods is only ${\sim}1\%-5\%$, even when the concentration of as-induced vacancies $(C_{0})$ by implantation is much higher than substitutional nitrogen [14, 29, 30]. One of the possible reasons for such low efficiency is that the rate of vacancy trapping into di-vacancies $\mathrm{(\textit{V}_{2})}$ or multi-vacancy complexes $\mathrm{(\textit{V}_\textit{n})}$ in such ion-damaged areas could be higher than trapping of single vacancy at substitutional nitrogen atom, which forms NV centre. At annealing temperatures in a range of 750 °C–950 °C, where vacancies become mobile [15, 21, 22, 29, 31], vacancies get trapped in $\mathrm{V_{2}}$ [32, 33] or $\mathrm{V_{n}}$ [17, 18]. Thus, the aggregation of vacancies should be considered as a competitive path for NV formation. Besides, if the formation energy of NVs is higher compared to that of $\mathrm{\textit{V}_{2}}$-related defects, the rate of NV formation at vacancy diffusion would be further suppressed [34].

To address this issue, we developed a model of vacancy diffusion based on probabilistic atomic jumps in crystal lattice, where the NV formation has been considered as competing path of vacancy trapping to that of $\mathrm{\textit{V}_{2}}$ and $\mathrm{\textit{V}_\textit{n}}$. One of the critical parameters of our model is a so-called capture cross-section, which can be related to the formation energy of each corresponding defect. Deriving an experimental estimates of the ratio of capture cross-sections for di-vacancies $(\sigma_{\textit{V}2})$ and NVs $(\sigma_\textrm{NV})$ is in focus of this study.

Vacancies have been induced in CVD and HPHT diamond crystals with low ppb nitrogen doping by helium implantation at different energies. The corresponding NV profiles have been generated upon annealing at temperatures of 950 °C and 1200 °C [14]. Precise evaluations of these NV profiles in depth have been obtained via 3D localization of individual NVs by a confocal microscope. In this way we evaluate the activation energy of vacancy diffusion in diamond. In turn it enables to evaluate the capture cross section for a di-vacancy.

Besides, electron exposure with MeV kinetic energies and in a wide range of doses $\mathrm{(10^{11}}$ to $\mathrm{10^{17}\, cm^{-2})}$ has been applied to evaluate the capture cross section ratio between NVs and $\mathrm{\textit{V}_{2}}$ for a homogeneous distribution of vacancies in the bulk. The results of this work can be interesting for further optimization of electron beam irradiation [31, 35] to create NVs while preserving their spin properties (bulk magnetometry).

The proposed model is based on probabilistic atomic jumps in a crystal lattice, which is a standard theoretical concept for describing atomic diffusion and recombination phenomena in solids [36]. As it has been shown by calculations, self-diffusion in diamond is dominated by vacancies due to the stiffness of C–C bonds in the lattice [37]. This is quite different from other group-IV semiconductors, where interstitials and direct-exchange mechanisms play a similar role. This situation simplifies our model significantly and facilitates the consideration of vacancy migration as a sequence of localized random jumps between neighboring sites in a 3D lattice. A 2D representation of the model is shown in figure 1(a). Here, a vacancy (white dot) jumps along four fixed directions (perpendicular to each other) exchanging position with carbon nearest neighbors (light-grey dots). During this process of sequential jumps of length a₀, this vacancy 'scans' a volume (an area shown in the 2D case) which then defines the probability of being captured by other defects, like by stationary impurities or by moving vacancies. This particular volume can be defined as $\sigma a_{0}$, where, $\sigma = \alpha a_{0}^{2}$ represents the capturing cross section of a defect and a₀ is defined as the interatomic spacing in the lattice.

Zoom In Zoom Out Reset image size

The parameter α is characteristic for every specific defect and can be explained in relation to its formation energy: in a hypothetical case, the energy of defect formation could be equal or even smaller compared to the energy barrier for diffusion of mono-vacancy in the lattice. According to this scenario, the placement of a vacancy in the defect neighboring position will suffice for trapping of the vacancy by the defect. Therefore, the capture cross-section of this particular defect would be equal to its geometrical capture cross section $(\alpha = 1)$. As shown in figure 1(a), it is given by a sum of all possible directions of the jumps into the neighboring cell. Each of the cell may include the targeted trapping centre. An example of vacancy jumps sequence is shown schematically in figure 1(a). The migrating vacancy is trapped at substitutional impurity atom (green dot) in the lattice, if both defects appear in the neighboring lattice sites. In this way a substitutional (S) impurity—vacancy complex (S-V) is formed, which for instance could be a NV centre in diamond.

It is reasonable to assume, that, for some defects the formation energy is higher than the activation energy of single vacancy jump. In this case, the mere placement of a vacancy in the defect neighboring position will not fulfill the necessary criteria for the vacancy to be trapped by the defect. The value of σ would therefore be below the geometric cross section and it would reflect a complex process of trapping which consists of several steps with different probabilities $(\alpha\lt1)$. Higher energy of defect formation leads to less probability for a vacancy trapping. Evaluation of such mechanism requires rigorous theoretical calculations which are beyond the scope of this study. Here, we refrain from replicating the underlying physics of a defect—vacancy complex formation phenomenon, but rather simulate it through a statistical approach. Thus, the value of a defect capture cross section is not crucial in our investigation. The relative value of a defect capture cross section compared to that of a vacancy is the key parameter of the model. We can assume that the physics governing the formation of NV centres is equivalent to that underlying the di-vacancy formation. Therefore, the $\sigma_\textrm{NV} / \sigma_{\textit{V}2}$ ratio evaluated in the model can be deemed valid not only for the model, but also for real-case scenarios. In this study the $\sigma_\textrm{NV} / \sigma_{\textit{V}2}$ ratio has been evaluated by applying the model results to experimental NV profiles.

The probability for a vacancy of being trapped during the diffusion (p) can be expressed as function of jumps (n) as

$$\begin{align} & p\left(n\right)\approx \dfrac{\sigma n a_{0}\left(1-\pi_{0}\right)}{\textrm{per}\: \textrm{unit}\: \textrm{Volume}} \end{align} \tag{ 1a }$$

$$\begin{align} \notag\\[-6pt] & n \longrightarrow n_\textrm{max}, \quad p\left(n\right) \longrightarrow 1 \end{align} \tag{ 1b }$$

where n is set to $n_\textrm{max}$ when the entire volume unit is completely 'scanned' by the jumping vacancy and the vacancy trapping probability approaches one. Here we take into account the probability $(\pi_{0})$ for a vacancy of returning to its previous lattice position after every jump. This probability depends on the lattice type, and for a vacancy in diamond is 0.44 [36].

In this work we describe vacancies trapping by competing defects for various cases of initial distributions of vacancies. In particular, we consider a homogeneous field of vacancies representing the case of high-energy electron exposure, and also vacancy clusters formed along the ion tracks during ion implantation. The reduction of free vacancies due to the formation of $\mathrm{\textit{V}_{2}}$, $\mathrm{\textit{V}_\textit{n}}$ and S-V will be addressed in the paper.

In order to preserve the probabilistic nature of the model, charge neutrality assumption is imposed on all vacancies and substitutional impurities. Moreover, an initial boundary of $C_{0} \gg S$ mimics the experimental conditions of low nitrogen doped substrates. Under this condition, the $\mathrm{\textit{V}_{2}}$ trapping is seemingly not affected by S-V formation. Hence, a hierarchy in the formation of the vacancies consuming defects has been introduced into the model, where the rate of vacancy aggregation is not affected by the formation of S-V defects. In our model, $\mathrm{\textit{V}_{2}}$ and $\mathrm{\textit{V}_\textit{n}}$ represent two extreme cases of defect formation in diamond in the given temperature range. $\mathrm{\textit{V}_{2}}$ model allows exact analytical description of the situation. Thus, it is used to check the validity of the assumption on the independence of $\mathrm{\textit{V}_{2}}$ and $\mathrm{\textit{V}_\textit{n}}$ from S-V trapping. This is made by comparing the analytical model to the results of Monte Carlo simulations (MC). Assuming a homogeneous distribution of initial 'as irradiated' vacancies, the formation of $\mathrm{\textit{V}_{2}}$ and $\mathrm{\textit{V}_\textit{n}}$ as function of jumps can be described as

$$\begin{align} dV_{2}\left(n\right) & = \sigma_{V_{2}}a_{0}\left(C_{0}-V_{2}\right)^{2}\left(1-\pi_{0}\right)dn \end{align} \tag{ 2a }$$

$$\begin{align} \notag\\[-6pt] dV_{n}\left(n\right) & = \sigma_{V_{n}}a_{0}C_{0}\left(C_{0}-V_{n}\right)\left(1-\pi_{0}\right)dn. \end{align} \tag{ 2b }$$

Here ${\textit{V}}_{2}$ and ${\textit{V}}_{\textit{n}}$ are the concentrations of vacancies which are trapped into di-vacancies and vacancy complexes respectively. In order to simplify the $\mathrm{\textit{V}_{\textit{n}}}$ model, free vacancies available to get trapped have been considered constant in this case. This condition might overestimate $\mathrm{\textit{V}_\textit{n}}$ formation with increasing the number of jumps n, in contrast to the case of $\mathrm{\textit{V}_{2}}$. On the other hand, the analytical approach of $\mathrm{\textit{V}_\textit{n}}$ defects allows an extension towards clustering of vacancies by ion implantation (see below). Analytical solutions of equations (2a ) and (2b ) can be evaluated as $V_{2}(n)$ and $V_{n}(n)$

$$\begin{align} V_{2}\left(n\right) & = \dfrac{\sigma_{V_{2}}a_{0}C_{0}^{2}\left(1-\pi_{0}\right)n}{1+\sigma_{V_{2}}a_{0}C_{0}\left(1-\pi_{0}\right)n} \end{align} \tag{ 3a }$$

$$\begin{align} \notag\\[-6pt] V_{n}\left(n\right) & = C_{0}\left\lbrace 1-exp\left[-\sigma_{V_{n}}a_{0}C_{0}\left(1-\pi_{0}\right)n\right]\right\rbrace . \end{align} \tag{ 3b }$$

The Einstein–Smoluchowski relation allows to relate the number of vacancy jumps in diamond lattice with the time $(t_\textrm{ann})$ and temperature $(T_\textrm{ann})$ of the annealing process [36] :

$$\begin{eqnarray} &&n_\textrm{jumps} = \dfrac{8t_\textrm{ann}D_{0}}{a_{0}^{2}\left(1-\pi_{0}\right)}\exp\left({-\dfrac{E_{a}}{k_{B}T_\textrm{ann}}}\right). \end{eqnarray} \tag{ 4 }$$

Here, k_B is the Boltzman's constant, E_a and D₀ represent the diffusion activation energy and the diffusion coefficient respectively. For diamond these parameters have largely been evaluated in literature as 1.7–2.1 eV [38, 39] and $3.69 \times 10^{-6}$ cm² s⁻¹ [32] respectively. In this work we made our own evaluation of E_a by applying the model to analyze experimental depth profiles of single NVs (see below).

Following the hierarchy of equation described above, the S-V defect formation can be expressed as

$$\begin{eqnarray} &&d\textrm{SV}\left(n\right) = \sigma_\textrm{SV}a_{0}\left(C_{0}-V_{x}\right)\left(S-\textrm{SV}\right)\left(1-\pi_{0}\right)dn, \end{eqnarray} \tag{ 5 }$$

where $\textrm{SV}$ and S define the concentration of S-V complexes and of substitutional defects respectively. V_x can refer to the concentration of di-vacancies $(x = 2)$ or of multi-vacancies $(x = n)$ complex as described in equation (3). The substitution of V₂ and V_n (evaluated through equations (3a ) and (3b ) in equation (5)) provides the following analytical solutions:

$$\begin{align} \textrm{SV}_{\left(V_{2}\right)}\left(n\right) & = S\left\lbrace 1-\left[1+\sigma_{V_{2}}a_{0}C_{0}\left(1-\pi_{0}\right)n\right] ^{-\dfrac{\sigma_\textrm{SV}}{\sigma_{V_{2}}}}\right\rbrace \end{align} \tag{ 6a }$$

$$\begin{align} \notag\\[-6pt] \textrm{SV}_{\left(V_{n}\right)}\left(n\right) & = S\left[1-\textrm{exp}\left(-\dfrac{\sigma_{SV}}{\sigma_{V_{n}}}\dfrac{V_{n}\left(n\right)}{C_{0}}\right)\right]. \end{align} \tag{ 6b }$$

Here, the σ parameters represent the corresponding capture cross sections for V₂, V_n and $\textrm{SV}$ respectively. $V_{n}(n)$ is given by equation (3b ) and describes the concentration of vacancies trapped into the $\mathrm{\textit{V}_{\textit{n}}}$-defects after n jumps.

In MC simulation vacancies depletion has been simulated by the formation of $\mathrm{\textit{V}_{2}}$, $\mathrm{\textit{V}_{\textit{n}}}$ and S-V defects. While $\mathrm{\textit{V}_{2}}$ and S-V can trap only one vacancy, no limitation on size has been set for $\mathrm{\textit{V}_{\textit{n}}}$ (C₀ is constant in equation (2b )). The calculation procedure by the MC simulation was terminated when all free vacancies got trapped. For a 3D cubic lattice, of which a 2D cut is shown in figure 1(a), we impose 6 possible directions of diffusion for a vacancy to move (cubic lattice for simplicity of calculation). The di-vacancy capture cross section corresponds to its geometrical size: $\sigma_{V2} = 6a_{0}^{2}$. The validity of this assumption has been verified through simulations that allow vacancies to move also along all diagonal directions, i.e. 26 degrees of freedom. These simulations involve varying the concentration of the initial free vacancies. Further details can be found in the supplementary material. Moreover, the assumption of neutral charge species establishes the same capture cross section for the S-V defect. For simplicity we set in the model the capture cross section of $\mathrm{\textit{V}_{\textit{n}}}$ as identical to $\mathrm{\textit{V}_{2}}$ and S-V.

In the first place vacancies trapping has been simulated for different concentrations of free vacancies allowing only $\mathrm{\textit{V}_{2}}$ formation. Figures 1(b) and (c) illustrate model (dashed curve) and simulation (solid curve) results of vacancies trapping due to $\mathrm{\textit{V}_{2}}$, $\mathrm{\textit{V}_{\textit{n}}}$ and S-V formation, under the initial condition $C_{0} = 10^{18}\, \mathrm{cm^{-3}}$ and $S = 10^{16}\, \mathrm{cm^{-3}}$. These concentrations have been simulated by including 5000 of initial vacancies $(C_{0})$ and 50 of P1-like impurities in a simulation volume, where the minimum distance between two neighboring atoms has been set as $a_\textrm{simulation} = 0.2\,\mathrm{nm}$ (more simulations show that the amount of jumps required to completely deplete the free vacancies decreases with the increment of the initial vacancies concentration. Information about the simulations are reported in the supplementary material). During the trapping process of multi-vacancy defects $\mathrm{V_{n}}$, it is clear that the amount of trapping sites reduces gradually with increasing the size of the vacancy chain. Some vacancies already trapped in the complex saturates their bonds and cannot acquire additional vacancies any more. This is the reason for the discrepancy in figure 1(c) for the case of $\mathrm{\textit{V}_{\textit{n}}}$ model vs. MC simulation. On the other hand, the agreement between model and MC simulation confirms our initial hypothesis that S-V formation can be considered as a negligible process that, when $C_{0} \gg S$, does not affect $\mathrm{\textit{V}_{2}}$ and $\mathrm{\textit{V}_{\textit{n}}}$ trapping.

However, in the real situation of implantation experiments, initial vacancies are formed in dense clusters due to collision cascades along the ion tracks at stopping. Such tracks are characterized by a high density of vacancies in the local environment. Depending on the mass and energy of the implanted atom, the local concentration of vacancies in clusters can range from ${\sim}10^{19} \: \mathrm{cm^{-3}}$ to $10^{21} \: \mathrm{cm^{-3}}$. Therefore, one would expect high loss of vacancies by trapping into vacancy complexes within such highly-dense but nanometer-size clusters by implantation. Thus, the number of vacancies available for NV activation will be reduced. According to the model equation (3b ), such complexes by local trapping are formed in a time scale much shorter compared to the diffusion and trapping in bulk of undamaged crystal associated with NV formation. This holds however only if the dose of implantation is still below a certain limit where neighboring implantation tracks do not coalesce.

For mathematical treatment of this problem it would be convenient to represent the total concentration of initial vacancies C₀ induced by implantation:

$$\begin{eqnarray} &&C_{0} = C_\textrm{free} + C_\textrm{free} \cdot K \end{eqnarray} \tag{ 7 }$$

where $C_\textrm{free}$ is a part of initial vacancies which 'escape' the border of the implantation-induced clusters and can contribute to activation of NVs by diffusion and trapping in the bulk of undamaged crystal. The term $C_\textrm{free} \cdot K$ describes here the other part of vacancies which are aggregated into complexes inside the volume of the initial cluster. The parameter K relates to the vacancy escape ratio which is $C_\textrm{free}/C_{0}$. As defined in this way it can be expressed as $1/(1+K)$.

To estimate the impact of such losses, we extend our MC and analytical models to the case of high density clusters of vacancies. Figure 1(d) shows the 2D schematic of the simulated structure. The vacancies are initially distributed inside a single cluster (blue circle) and can aggregate within the cluster or leave its borders (white dots). No limitation was imposed on the amount of vacancies in the aggregate. In figure 1(e) simulation ends when there are no remaining free vacancies within the initial cluster's volume. Simulations were performed to investigate the vacancy escape ratio from clusters with varying their vacancy concentration. This has been realized by modifying the cluster size $(d_{0})$ for clusters initially consisting of ten (purple dots) and fifty (green dots) vacancies. The size of the cluster was varied in order to correlate with the density of vacancies within the stopping range by He, C and N implantation (10¹⁹–$10^{21} \: \mathrm{cm^{-3}}$ per cluster in average). The results of MC simulation in figure 1(e) were fitted by equation (2b ) in order to estimate the number of jumps required to fulfill the conditions of vacancy escape in each particular case. Here we found a kind of a linear proportionality between the average diffusion path of vacancies to escape from a cluster, i.e. the minimum number of jumps, and the cluster size: $\langle d_\textrm{esc}\rangle \sim kd_{0}$. A fit by the model from equation (2b ) is shown in figure 1(e) by the dashed curves, yielding a proportionality constant k ≈ 0.33. We incorporate these results on the escape ratio into our model to check this effect on the yields of S-V defects.

In the next step we perform simulations on S-V formation considering the clusters of initial vacancies. The results of modelling were fitted by analytical model of multi-vacancy defects $(\mathrm{\textit{V}_{\textit{n}}})$, modified in accordance with the discussion above. Similar to the case of homogeneous distribution by equation (2b ), the rate of $\mathrm{V_{n}}$ formation outside the implantation clusters can be described as:

$$\begin{eqnarray} &&dV_{n, <\textrm{cluster}>}\left(n\right) = \sigma_{V_{n}}a_{0}C_\textrm{free}\left(1+K\right)\left(C_\textrm{free}-V_{n}\right)\left(1-\pi_{0}\right)dn. \end{eqnarray} \tag{ 8 }$$

The index $\lt \textrm{cluster}\gt$ here is to distinguish the cluster-related model from the case of homogeneous distribution by equation (3b ). The parameter K describes trapping inside the implantation-related clusters according to equation (7). For each particular implantation, the vacancy escape ratio is given by the simulation results in figure 1(e). Then the corresponding parameter K and $C_\textrm{free}$ can be calculated. For instance, for a typical escape ratio of 0.2 the value of K is 4.0. In equation (8) we set that free vacancies $C_\textrm{free}$ escaped from the aggregation inside clusters can diffuse in the lattice and participate in the formation of $\mathrm{\textit{V}_{2}}$, $\mathrm{\textit{V}_{\textit{n}}}$ and S-V defects. This is described by the term $C_\textrm{free}-V_{n}$. On the other hand we impose no limitation on the amount of vacancies in the aggregate. Therefore, these free vacancies can also diffuse into the volume of neighboring clusters and get trapped there at the immobile vacancy complexes. Then the total amount of trapping centres in equation (8) is given by the term $C_\textrm{free}(1 + K)$, which is the total concentration of the initial, i.e. as implanted vacancies C₀ by equation (7).

Under these assumptions, integration of equation (8) yields the following analytical solution:

$$\begin{eqnarray} &&V_{n, <\textrm{cluster}>}\left(n\right) = C_\textrm{free}\left\lbrace 1-\textrm{exp}\left[-\sigma_{V_{n}}a_{0}C_\textrm{free}\left(1+K\right)\left(1-\pi_{0}\right)n\right]\right\rbrace. \end{eqnarray} \tag{ 9 }$$

Note, that equation (9) describes vacancy trapping into complexes outside the implantation-related clusters. Only this part is relevant for the S-V formation. The formation of S-V defects in the cluster model can be expressed then following the same hierarchy of equations as in the case of homogeneous distribution:

$$\begin{eqnarray} &&\textrm{SV}_{<\textrm{cluster}>}\left(n\right) = S\left[1-exp\left(-\dfrac{\sigma_\textrm{SV}}{\sigma_{V_{n}}}\dfrac{V_{n,<\textrm{cluster}>}\left(n\right)}{C_\textrm{free}\left(1+K\right)}\right)\right]. \end{eqnarray} \tag{ 10 }$$

where $V_{n,\lt \textrm{cluster}\gt}(n)$ is expressed by equation (9), compare with equation (6b ).

Figure 1(f) shows the dependence of the S-V formation efficiency on the vacancies cluster escape ratio (blue curve). Here, simulation (dots) and model (dashed curve) results are characterized by $C_{0} = 10^{17} \,\mathrm{cm^{-3}}$ and $S = 10^{16}\, \mathrm{cm^{-3}}$, with vacancies initially distributed in clusters where $C_\textrm{cluster} \gg C_{0}$. As expected, the presence of vacancies clusters related to ion implantation imposes an additional reduction in the S-V formation. Therefore, it is evident that opting for the implantation of low mass atoms is a better choice.

According to the model's prediction, if the NV centre is the S-V defect, the efficiency of the P1 to NV conversion is still considerably lower than what has been observed in experiments $({\sim}1\%$–5%). This discrepancy could be related to different capture cross sections per $\mathrm{\textit{V}_{2}}$ and $\mathrm{\textit{V}_{\textit{n}}}$ formation.

In the experimental section we apply our model to the experimental data on electron and ion irradiation of diamond in order to get estimates on the capture cross-sections for NV and $\mathrm{\textit{V}_{2}}$ formation.

3.1. Sample preparation

3.2. 3D optical localization of NV centres for depth profiles evaluation

3.3. A near-surface depletion layer of NV centres

NV depth scans of electron-exposed diamonds at low-dose limit revealed the presence of a NV sub-surface depletion region (figure 2). An increase in the exposure dose leads to a gradual reduction of the depletion region. For this particular sample shown in figure 2 the depletion layer disappears at a dose of e-exposure (10 MeV) of $1 \times 10^{16}$ el cm⁻². Tentatively we relate it to sub-surface damage caused by sample polishing. Similar effects have been reported in literature where the appearance or the increase of defect bands in the cathodoluminescence spectra have been removed by plasma etching [42, 43]. Such structural defects, like dislocation, could trap single vacancies more effectively than nitrogen donors, thus leading to the NV-depletion region. An increment of the exposure dose above a certain limit could saturate all competing trapping centres.

Zoom In Zoom Out Reset image size

A schematic representation of the trapping effect is shown in figure 2(d). To minimize a potential depletion effect by polishing, samples used in this study for NV depth profiling were subjected prior to implantation to a RIE/ICP plasma in Ar:SF6 atmosphere in order to etch approximately $3\,\mu \mathrm{m}$ in depth. Chemical etching in hot piranha solution has been applied after Ar:SF6 etching to restore the O-terminations at the surface.

The PL scan of the substrate after e-irradiation at high dose in figure 2(c) reveals a layered pattern of NVs in depth. Such distribution profiles are typically observed in nitrogen assisted CVD growth [44]. This creates a stepped surface which provides different crystallographic planes that have different uptake in N and causes striations in NV luminescence through the diamond.

4.1. NV profiles by helium-implantation and annealing

In this work we study the evolution of the NV depth profiles by He-implantation in the temperature range of 950 $^{\circ}\mathrm{C}$–1200 $^{\circ}\mathrm{C}$. Thermal annealing not only induces vacancies diffusion, but, leads to different type of defects in irradiated diamond depending on the applied temperature range [14]. A reduction of single vacancies up to $30\%$ is observed after annealing at 450 $^{\circ}\mathrm{C}$ [32, 33, 45, 46]. Single vacancies, immobile at this temperature, recombine with mobile interstitial defects. We account for this effect when estimating the concentration of free vacancies in our model. The main trapping part above 600 $^{\circ}\mathrm{C}$ is the vacancy aggregation in di-vacancies complexes [32, 33]. Efficient NV formation occurs at temperature of 800 $^{\circ}\mathrm{C}$ or above. A further increment of the temperature up to 1000 $^{\circ}\mathrm{C}$ leads a di-vacancies disappearance to the benefit of multi-vacancies aggregates formation [17, 18]. At such conditions, the NV formation can be considered as a competitive trapping of vacancies between nitrogen atoms and vacancies complexes.

The experimental NV depth profiles by helium irradiation in a 100-oriented CVD diamond $(\mathrm{P1} \sim \mathrm{1ppb})$ are shown in figures 3 and 4. For the experiments in figure 3, the initial condition $C_{0} \gg \mathrm{P1}$ has been ensured by a helium dose of $5 \times 10^{11}\,\mathrm{cm^{-2}}$, at 7 °-off angle. This produces approx. $1.5 \times 10^{18}\,\mathrm{cm^{-3}}$ of vacancies at the ion stopping range according to Crystal Trim (C-TRIM) simulation with a displacement energy of 35 eV. The areal densities of NVs in figures 3 and 4 have been converted to a concentration considering the scanning area of $10\,\mu \mathrm{m} \: \times \: 10\,\mu \mathrm{m}$ and the depth steps between 2 consecutive scans.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The NV profile after annealing at 950 $^{\circ}\mathrm{C}$ for 2 h shown in figure 3(a) (blue histogram) follows the vacancies distribution by C-TRIM for the same experimental conditions. Here, it is necessary to mention, that NV depth by He irradiation and 950 $^{\circ}\mathrm{C}$ annealing is limited mostly by local trapping of vacancies within the initial ion-implanted region (see also figure 6 discussed below).

This conclusion is also confirmed by the statistics of $\mathrm{T_{2}}$ (Hahn-echo) measurement taken at the corresponding NV profiles. The Hahn-echo sequence was applied at an external magnetic field of 120G. As shown by the blue dots in figure 3(b), approx. $30 \%$ of the total analyzed NVs show a $\mathrm{T_{2}}\lt20\, \mu \mathrm{s}$ and a noticeable component of $\mathrm{NV^{0}}$ in the PL spectra. On the other hand, a spin coherence time up to ${\sim}340\, \mu \mathrm{s}$ can be observed for NVs distributed most probably deeper in the bulk indicating a reduction in the lattice damage. The similarity between the NV profile and the helium stopping range, along with the low $\mathrm{T_{2}}$ suggests that trapping of vacancies due to aggregation in the vicinity of radiation defects (i.e. within the helium stopping range) occurs at a faster rate than the projected macroscopic diffusion. This result is in accordance with the model which predicts a high density of di-vacancies and multi-vacancies complexes within the stopping range of helium atoms.

A subsequent annealing of sample A at 1200 $^{\circ}\mathrm{C}$ for 1 h generates the NV profile shown in red histogram in figures 3(a) and 4. Note, that the two scans after 950 $^{\circ}\mathrm{C}$ and 1200 $^{\circ}\mathrm{C}$ annealing steps were made at approximately the same area of the sample. Similar NV profiles at the same conditions for another sample (sample B) are shown in figure 4(a). Here the NV depth profiles collected from two different regions are characterized by a depth ${\gt}1\, \mu \mathrm{m}$, which is a demonstration of the long range diffusion of vacancies from the initial irradiation damaged area.

As a result of higher temperature annealing, three effects are observed: (i) the vanishing of NVs with low $\mathrm{T_{2}}$ (see figure 3(b)), (ii) the deeper profile of NVs into the bulk, and (iii) the inhibition of pre-existing NVs in the region damaged by implantation (two histograms in figure 3(a)).

The disappearance of NVs with $\mathrm{T_{2}}\lt20\,\mu \mathrm{s}$ suggests a reduction of lattice damages within the He stopping range, which has been achieved through elimination of residual implantation related paramagnetic defects and vacancy chains dissociation [14]. Additionally, an improvement of the spin coherence time up to $500\,\mu \mathrm{s}$ is found for NVs created in the deep tail (depth ${\gt}1\,\mu \mathrm{m}$). This indicates that no competing formation of $\mathrm{\textit{V}_{\textit{n}}}$ defects takes place (at least in a range of detectable concentration, like in the stopping range at 950 $^{\circ}\mathrm{C}$, where $\mathrm{\textit{V}_{2}/\textit{V}_{\textit{n}}}$ prevail). The initial irradiation damaged area, therefore, operates as a vacancies source. Here, indeed, single vacancies are released from chains and diffuse into the bulk to a depth ${\gt}1\,\mu \mathrm{m}$. The longer vacancies range diffusion is explainable with an increase of the amount of jumps per single vacancy compared with annealing at 950 $^{\circ}\mathrm{C}$ for 2 h as described by equation (4).

The increase of the NVs depth distribution and of the spin coherence time for annealing at 1200 $^{\circ}\mathrm{C}$ suggests that vacancy diffusion occurs much faster than vacancy aggregation in the deep tail. Hence, in this case the formation of vacancy complexes can be considered negligible compared to that of NVs. Based on these arguments, we modified our model in order to extract the activation energy of vacancy diffusion by fitting the deep tails of the NV distribution. We consider the case of a remote source of vacancies, which is located within the stopping range of He atoms. In a region away from this source, the concentration of NVs at each point of depth will be given by the product:

$$\begin{eqnarray} &&\textrm{NV} \left(x,t\right) \sim f\left\lbrace \sigma_\textrm{NV} \left[P1\right] \int C\left(x,t\right)dt\right\rbrace, \end{eqnarray} \tag{ 11 }$$

where [P1] is the local concentration of substitutional nitrogen donors taken as constant through all depth. Here, the vacancy trapping on nitrogen donors is assumed to be not limited by multivacancy defects, as in equations (6a ) and (10). The $C(x,t)$ function depends on the chosen model of vacancy source. We consider two opposite cases: instantaneous or constant [36]. An instantaneous source is characterized by a 'flash-like' emission of vacancies, which is much shorter than the duration of the annealing process. This emission can be described as

$$\begin{eqnarray} &&C\left(x,t\right) \propto \textrm{exp} \left[\frac{-x^{2}}{4D_{0}\textrm{exp}\left(-E_{a}/k_{B}T_\textrm{ann}\right)} \right]. \end{eqnarray} \tag{ 12 }$$

A constant source, instead, is characterized by a vacancies emission comparable with the duration of the annealing process:

$$\begin{eqnarray} &&C\left(x,t\right) \propto erfc \left(\frac{x}{2\sqrt{D_{0}\textrm{exp}\left(-E_{a}/k_{B}T_\textrm{ann}\right)t}} \right). \end{eqnarray} \tag{ 13 }$$

The results of the model are shown by dashed lines in figure 4(a). Both approaches yield a vacancy diffusion activation energy of about 1.7 eV, which is quite consistent with earlier estimates [38]. This value is close to the theoretical estimates of energy barrier for di-vacancies formation in diamond [47] (${\sim}1.8 \mathrm{eV}$ of the lowest path in the calculation). Therefore, as previously discussed, we can consider the capture cross section for the di-vacancy defect to be close to the geometrical one, $\sigma_{V2} = 4a_{0}^{2}$. The factor 4 stands here for the number of bond directions in $\text{sp}_{3}$ lattice. It is also worth to note, that the use of the equation (12) or (13) to fit the E_a directly from NV profiles might not be correct: such functions describe, indeed, only the diffusion of vacancies. The resulting NV profiles in our study, instead, is a product of two separate processes: vacancy diffusion and vacancy trapping on P1 at each point of depth.

As mentioned before, the additional annealing at 1200 $^{\circ}\mathrm{C}$ not only generates long vacancies diffusion, but, also induces a reduction in the concentration of the NVs already formed in the area damaged by ion implantation. Similar NVs quenching phenomenon has been widely reported in literature [14, 15, 48, 49] as a consequence of annealing at temperature higher than 1000 $^{\circ}\mathrm{C}$. One possible reason for NV quenching could be the formation of NVH defect. To prevent this defect from depleting NVs, samples A and B were pre-annealed at 950 $^{\circ}\mathrm{C}$ for 2 h before ion implantation. This pre-annealing step aimed to remove any residual atomic hydrogen from the lattice after CVD growth process. On the other hand, the increment of the spin coherence time measured in the ion damaged region suggests the absence of states in the energy bandgap that could electrically compensate NVs. Indeed, defects which induce electronic states in the bandgap for charge trapping would also lead to a reduction of electron spin coherence of the NV centres in the vicinity.

A possibility of vacancies surface recombination has been excluded through NV depth distributions generated by low energy ion implantation (figure 4(b)). Electronic grade $(\mathrm{P1}\sim \mathrm{1ppb})$ CVD diamond (sample C) has been Helium implanted with a dose of $2 \times 10^{11}$ cm⁻² and $5 \times 10^{11}$  cm⁻² respectively at 2.5 eV and 5 eV. A micro-beam ion implanter with coaxial optics for sample alignment has been used, with a nominal-off angle of approx. 7 °. The low implantation energy provides, in accordance to C-TRIM simulation, a vacancies distribution restricted to 10–20 nm below the surface. Figure 4(b) displays the NV depth profiles obtained from sample C following the same annealing procedure as samples A and B. These profiles reveal a formation of NVs that extend up to $5\,\mu \mathrm{m}$ into the bulk, providing further evidence that vacancies do not exhibit a preference for surface recombination. Furthermore, the deep tails of the NV profiles can be accurately modelled using an activation energy of 1.7 eV.

Thus, the reason behind the reduction of NV density could be related to single vacancies which are released at 1200 $^{\circ}\mathrm{C}$ from previously formed $\mathrm{V_{n}}$ defects, and which diffuse deep in the bulk thus forming single NVs visible in the depth tails in figures 3 and 4. These vacancies bind with NVs to form the $\mathrm{N\text{-}V\text{-}V}$ defect [50, 51]. Alternatively, near-surface traps, might bind electrons from NVs and render them unobservable by optical means. The presence of traps combined with different capturing cross section per defect could explain the low P1 to NV conversion efficiency associated to ion implantation technique.

4.2. NVs in bulk by e-exposure: capture cross-section of vacancy trapping

Once the diffusion activation energy is verified, the ratio of $\sigma_{NV}/\sigma_{V2}$, previously discussed, can be deduced through the experiments on electron exposure. Deep in the bulk NVs are not influenced by any surface-related effects and are far from the depletion region (see figure 2). For this purpose we irradiated several 100 oriented CVD diamond substrates (1ppb) with 10 MeV at different doses.

Samples were irradiated with a 10 MeV linear electron accelerator (MB10-30MP; Mevex Corp.), where the dose was determined by means of a graphite calorimeter. The electron fluence then was re-calculated taking into account the mass collision stopping power in the corresponding energy range.

To estimate the density of as-irradiated vacancies, we used data from reference [52] for vacancy conversion of 3.42 vacancies per electron cm⁻¹ for 10 MeV electrons in diamond. This has been further used to evaluate the initial amount of vacancies produced by a particular dose, and which was considered as reduced by $30\%$ to consider the effect of vacancy-interstitial recombination at lower temperatures of about 400 $^{\circ}\mathrm{C}$ [32, 33, 45, 46]. We estimated the NV concentration deep in the bulk below the NV depletion region by 3D PL mapping technique, as for the NV profiles discussed above. For elevated electron doses, the NV concentration has been evaluated by comparing the photoluminescence intensity collected from such maps with the one of a single NV. The measured NV concentrations have been fit by the model considering only the NV and $\mathrm{V_{n}}$ formation. Equation (6b ) has been used for the case $C_{0} \gg \mathrm{P1}$. In the case of low dose exposure, $C_{0} \ll \mathrm{P1}$, instead, the effect of $\mathrm{V_{n}}$ formation can be considered negligible compared to that of NV ($V_{x} = 0$ in equation (5)). Under this assumption the rate of NV formation can be described as

$$\begin{eqnarray} &&\textrm{SV}\left(n\right) = S\left( 1-\textrm{exp}\left[-\sigma_\textrm{SV}a_{0}C_{0}\left(1-\pi_{0}\right)n\right]\right). \end{eqnarray} \tag{ 14 }$$

As shown in figure 5 by the dashed curves, the fit provides a NV capture cross section of $0.1\sigma_{V_{n}}\lt \sigma_\textrm{NV} \lt0.5\sigma_{V_{n}}$, where the vacancy chain capture cross section has been set as the geometrical one, $\sigma_{V_{n}} = 4a_{0}^{2}$. In addition to the data on low-ppm CVD diamond, figure 5 includes also the data on electron exposure of CVD diamond crystals with nitrogen doping in the range of 1 ppm. These samples were produced by CVD in a microwave-plasma-assisted synthesis with 100 major faces containing neutral substitutional nitrogen of approx. 1 ppm. The nitrogen concentration was determined by ultraviolet-visible absorption measurements at 77 K [53]. Different irradiation doses applied to each sample within a range of $10^{17}\,\mathrm{cm^{-2}}$ showed NV concentrations ranging from 180 to 430 ppb.

Zoom In Zoom Out Reset image size

The results of the modelling establish that vacancies aggregation in multi-vacancy complexes occurs with higher probability than the NV formation. It can be noted here that the conclusion on the low capture cross-section of NV trapping agrees also with the results of theoretical calculations on defect formation energies in diamond [34]. Thus, one would require unavoidably a high dose of electron exposure to activate NVs in diamond bulk by irradiation. Only a small portion of vacancies will be trapped on P1, while the rest will form $\mathrm{\textit{V}_{\textit{n}}}$.

4.3. NV by helium implantation. Limits due to near-surface traps

The obtained data on the vacancies diffusion activation energy and the capture cross section allows us to describe the NV formation within the near surface region. As already discussed, an example of such analysis is shown in figure 6, where a CVD grown electronic grade $(\mathrm{P1} \sim \mathrm{1ppb})$ diamond (sample D) has been implanted with Helium at 340 keV, with a dose of $1 \times 10^{11}$ cm⁻² and subsequently annealed at 950 $^{\circ}\mathrm{C}$ for 2 h. The dose was chosen in order to have a minimum radiation damage but still affordable NV statistics in such low nitrogen doped substrate. The NV depth profiles for different offset implantation angles, (nominally 0 ° and 7 °) are shown by histograms in figure 6. According to the model, even assuming a partial loss of vacancies by trapping within the implantation-related clusters and the lowest $\sigma_\textrm{NV}-to-\sigma_{V2}$ ratio of 0.1, the P1 to NV conversion efficiency should be ${\sim}10\%$. The reason behind the relatively low NV yield might be related to addition effect of near-surface traps, responsible for subsurface NV depletion region of electron exposed diamond crystals at low doses (see figure 2). To estimate the traps concentration, we introduce them in our model, assuming their capture cross section identical to that of $\mathrm{\textit{V}_{\textit{n}}}$ $(\sigma_\textrm{traps} = \sigma_{V2} = \sigma_{\textit{V}_{\textit{n}}} = 4a_{0}^{2})$. The introduction of traps modifies the NV formation equation such as

$$\begin{eqnarray} &&\textrm{SV}_{<\textrm{cluster}>}\left(n\right) = S\left[1-\textrm{exp}\left(-\dfrac{\sigma_\textrm{SV}}{\sigma_{V_{n}}}\dfrac{V_{n,<\textrm{cluster}>}\left(n\right)}{\left(1+K\right)\left(C_\textrm{free}+C_\textrm{traps}\right)}\right)\right] \end{eqnarray} \tag{ 15 }$$

where $V_{n, \lt \textrm{cluster}\gt}(n)$ represents the concentration of vacancies in trapped into vacancy complexes after n number of jumps. Fitting the experimental data yields the concentration of about $2 \times 10^{17}$ cm⁻³ when $\sigma_\textrm{NV} = 0.1\sigma_{V_{n}}$. The results on the traps modelling are summarized in figure 5 as well. It shows the experimental data on 340 keV He irradiation in nominally 1-ppb CVD diamond with as-polished surface, and also the data on 4 keV Helium irradiation taken from our earlier work in [29]. Here we include also the data on the complimentary measurement on HPHT diamond crystal with approx. 3 ppb nitrogen doping exposed to 4 keV Helium at the similar range of doses as in [29]. The sample was grown by the temperature gradient method at 5.5 GPa and ${\sim}1350\,^{\circ}\mathrm{C}$ by using Co-Fe solvent with Ti added as a nitrogen getter. The (100) plate was obtained by laser-cutting and polishing with diamond fine powders. For quantifying the low concentration of nitrogen donors, rapid-passage technique of continuous wave electron paramagnetic resonance (cw-EPR) was employed [54].

Zoom In Zoom Out Reset image size

The C₀ in this graph is related to the maximum concentration of vacancies by C-TRIM simulation. The final value shown on the graph had been reduced by $30\%$ to account for vacancy-interstitial recombination as for the case of e-exposure as discussed above. Regarding the data points on NV-to-P1 ratio for 4 keV experiments, we used the peak concentration of NVs as evaluated from the NV-depth profiling at some irradiation doses in [29]. In summary, the dash-dotted curves 'Traps' in figure 5 shows the fitted NV conversion ratio for a traps concentration of approx. $2 \times 10^{17}$ cm⁻³ when $\sigma_{NV} = 0.1\sigma_{V_{n}}$. The derived value of traps concentration is in a similar range as that estimated from the dose limits shown in figure 2. It is possible to observe that the trap concentration fits approximately all data points related to CVD and HPHT diamonds. The origin of these traps could be related to an initial surface polishing procedure for both CVD and HPHT substrates. In these regards, the method of NV depth profiling in PL scans might help to identify their origin and to find a way to handle it in the future.

In this paper we investigated some aspects of vacancies diffusion in diamond lattice with focus on the efficiency of NV formation though vacancy trapping at the nitrogen donor atoms (P1 centres) upon thermal diffusion. To study the limits of the P1-to-NV conversion by this technique we developed an atomistic model of vacancy trapping in a crystal lattice. This model is based on probabilistic atomic jumps, considering competitive trapping of mono-vacancies between di-vacancies $(\mathrm{\textit{V}_{2}})$, or multi-vacancy complexes $(\mathrm{\textit{V}_{\textit{n}}})$ and targeted NVs. In the model, we defined the probability of a defect to trap a vacancy as the defect capture cross section. The validity of this approach was verified by Monte-Carlo simulation of various trapping scenarios in a 3D lattice.

The results of modelling were applied to experimental NV profiles generated in CVD and HPHT diamond substrates upon electron exposure, or helium implantation and subsequent thermal annealing. In particular, several electronic-grade CVD diamond substrates were implanted with helium ions at 340 keV and also 4 keV energies and subsequently annealed at 950 $^{\circ}\mathrm{C}$ and 1200 $^{\circ}{C}$ temperatures. The resulted NV depth profiles were analyzed using 3D localization of individual NVs in a confocal microscope and sampling their spin coherence properties $(\mathrm{T_{2}})$ to identify the local environment of trapping centres. In addition, we measured the NV density in the bulk of CVD diamonds with nitrogen doping at low-ppb and low-ppm levels, subjected to MeV electron exposures in a wide range of doses.

Applying the model to the experimental data we estimated the activation energy of vacancy diffusion in diamond about 1.7 eV, which is close to earlier estimates in the literature. The ratio of defect capture cross sections related to di-vacancy $(\sigma_{V2})$ and NV $(\sigma_{\textrm{NV}})$ and defined as $\sigma_\textrm{NV}/\sigma_{V2}$ was about 0.1–0.5. This result demonstrates the tendency of vacancies to aggregate together rather than being trapped by nitrogen, defining this way a limit of the classical approach in the NV formation by ion implantation or electron irradiation followed by annealing.

In addition, details on NV depth profiling by helium irradiation allowed to get an insight into dynamic of vacancy trapping in diamond. It has been confirmed that vacancy aggregation within ion stopping range at about 950 $^{\circ}\mathrm{C}$ annealing dominates over macroscopic diffusion due to high density of prime radiation damages. During the subsequent annealing at 1200 $^{\circ}\mathrm{C}$ these aggregates can be regarded as a source of vacancies which can diffuse to a depth of several µm and activate NV centres by trapping. Due to the intact lattice environment, the spin coherence time of such deep NVs ranges up to $500\,\mu \mathrm{s}$ in the used CVD substrates.

Finally, we observed two subsidiary effects: (i) the density of NVs within the ion stopping range by helium implantation was noticeably reduced upon 1200 $^{\circ}\mathrm{C}$ annealing. At the same time, the NVs with low $\mathrm{T_{2}}$ below approx. $10\,\mu \mathrm{s}$ disappear. This effect can be related to formation of $\mathrm{N\text{-}V\text{-}V}$ defects formed by release of single vacancies from the initial vacancy aggregates located at the nanometer proximity to the 'low-$\mathrm{T_{2}}$' NVs. Further investigation of this effect could be beneficial in understanding the origin of a low yield of NV centres by low-energy nitrogen implantation.

(ii) Another effect was related a subsurface NV depletion region observed in the electron irradiated samples at low dose of irradiation. Since the recombination of vacancies at the surface as the main path of extinction was experimentally excluded, we proposed alternative traps in the model to explain the presence of such depletion region. The origin of traps could be related to pre irradiation polishing procedure, however reliably identifying the cause requires further investigation.

The presented model, by virtue of its probabilistic approach, could be used in the future as a tool for engineering the production of other defects in diamond.

This work has been supported by the DFG via the FOR2724 and the GRK2624, by the BMBF via Qsens and by EU via the AMADEUS Project.

All data that support the findings of this study are included within the article (and any supplementary files).