Splitting of photoluminescent emission from nitrogen–vacancy centers in diamond induced by ion-damage-induced stress

We report a systematic investigation on the spectral splitting of negatively charged, nitrogen–vacancy (NV−) photoluminescent emission in single-crystal diamond induced by strain engineering. The stress fields arise from MeV ion-induced conversion of diamond to amorphous and graphitic material in regions proximal to the centers of interest. In low-nitrogen sectors of a high-pressure–high-temperature diamond, clearly distinguishable spectral components in the NV− emission develop over a range of ∼4.8 THz corresponding to distinct alignment of sub-ensembles which were mapped with micron spatial resolution. This method provides opportunities for the creation and selection of aligned NV− centers for ensemble quantum information protocols.

centers [22,23]. Since the 1960s, experiments have been conducted to study the effect of uniaxial mechanical stress along the 100 , 110 and 111 directions on the spectral properties of the zero-phonon-line (ZPL) absorption and photoluminescence (PL) emission of characteristic defect-related centers in diamond, such as for example the neutral vacancy GR1 transition at λ = 741 nm [24][25][26][27][28], the neutral nitrogen-vacancy (NV 0 ) transition at λ = 575 nm [29,30], the H3 line at λ = 503 nm [31,32] and the H4 vibronic band at λ = 496 nm [33,34]. In particular, absorption measurements of the ZPL NV − transition under uniaxial stress were performed allowing the attribution of this band to a transition between the A 1 (ground) and E (excited) electronic states of a trigonal center [12] and more recently to explore the 1046 nm transition [35]. More recent work has been focused on the investigation of the effect of local strain on the splitting of the excited levels structure of the NV centers with both experimental measurements and first-principles modeling of the excited levels structure of the center [36][37][38].
With a thorough understanding of the effect of strain on the splitting and broadening of the ZPL emission, the opportunity to use these effects to monitor internal stresses in artificial diamonds now arises. In flame-grown polycrystalline diamond films produced by chemical vapor deposition (CVD), the splitting of the ZPLs of the 533 nm and the NV 0 emissions was directly measured and mapped across the cubo-octahedral crystallites using cathodoluminescence (CL) [39]. The 533 nm emission was studied in greater detail, while the NV − emission was not acquired since the center is not CL-active [40]. The effect of stress on PL emissions from NV 0 (λ = 575 nm) and Si (λ = 738 nm) was also correlated with electron spin resonance measurements from the substitutional nitrogen center (P1) in both undoped and nitrogen-doped CVD diamond films [41]. Similarly, the linewidth of the GR1 neutral-vacancy-related emission at 741 nm was employed to qualitatively monitor the strain in a range of brown-colored type IIa diamonds before and after thermal annealing [42]. Stress-related ZPL shifts of the NV 0 emission were also recently measured in diamond nanocrystals using CL micro-mapping [43].
Here, we report the direct observation of the splitting of NV − PL emission engineered by the mechanical stress induced by MeV ion implantation. The implantation of energetic ions into diamond induces damage and amorphization, with associated swelling [44]. Here we show that the intense stress fields generated by such distortions induce a splitting in the NV − PL emission that can be mapped with micrometric spatial resolution in regions extending beyond the implanted regions. In the present experiment, PL mapping was performed at cryogenic temperatures (T = 80 K) on type Ib diamonds which had been implanted with 2 MeV He + ions, both before (i.e. as-implanted) and after thermal annealing. This allowed the direct observation of stress-induced splitting of the NV − emission at variable distances from the implanted region, which was also correlated with both finite element method (FEM) simulations of the mechanical stress and Raman experimental data.
Although there are clear applications for our results in protocols using the ground-state properties of the NV − center, for example magnetometry [9] and electrometry [45], we have not directly measured the change in the ground-state properties as a function of the implantationinduced strain. As the Stark shift for the ground state is one order of magnitude smaller than that of the optical transition, the perturbations to the ground-state Hamiltonian are likely to be large, and would need to be taken into account in any practical device design.

Experimental
The sample under examination was a 3.5 × 3.5 × 1.5 mm 3 type Ib (10 ppm [N] 100 ppm) synthetic high-pressure-high-temperature (HPHT) single-crystal diamond (Sumitomo). The sample was cut along the 100 crystal direction and optically polished on the two opposite large faces. Artificial diamonds produced by this technique are structured in macro-sectors that develop during the HPHT growth from a single seed; depending on their orientation, they are characterized by different impurity concentrations (N, Ni, etc) and consequently they display distinct luminescence features [46]. As shown in the CL map of figure 1, the sectors are relatively large and structured along the crystallographic directions. In comparison with the 'main' sector (A), the 'strip' sectors extending near the sample edges (B) are characterized by lower nitrogen concentration, while the sectors near the sample corners (C) contain high concentrations of impurities (N, Ni, Fe, etc), as confirmed by PL mapping [47,48].
It is worth noting that although an HPHT sample was chosen with the aim of providing a strong signal in micro-PL measurements due to the relatively high concentration of native nitrogen, particularly significant results were obtained from low-nitrogen sectors (see below), thus indicating that stress-engineered splitting in single NV centers has a significant potential in high-quality CVD diamonds.
The sample was implanted with 2 MeV He + ions at the MP2 microbeam line of the 5 U NEC Pelletron accelerator at the University of Melbourne. A region of 280 × 125 µm 2 extending across the 'main' and 'strip' sectors was irradiated in the region highlighted by the red square in figure 1, as shown in figure 2. An implantation fluence of 1 × 10 17 cm −2 was uniformly delivered by raster scanning the ion microbeam; the ion current was ∼2 nA and the sample was Optical transmission microscopy image of the implanted region (highlighted in the red oval) extending across the 'main' sector (bottom) and the low-nitrogen 'strip' (highlighted in blue). The bottom-left inset shows a zoom image of the implanted area, where it can be appreciated that the masked side of the implantation (right side) is sharper than the unmasked side (left side). The top-right inset shows the schematic diagram of the implanted region and of the 'strip' sector which is employed in the following figures to highlight mapped regions. The blue arrow highlights the direction of the PL scan to locate the 'strip' sector reported (see figure 10).
implanted at room temperature. Using the SRIM 2008.04 Monte Carlo simulation code [49], we estimated the linear damage profile, expressed as the number of vacancies formed per incoming ion at a given depth. The simulation was performed by setting the atomic displacement energy to 50 eV [50] and, as shown in figure 3, the structural damage is mainly concentrated at the end of the ion range in the material (∼3.5 µm). A crude estimation of the damage density in the heavily damaged layer can be obtained by multiplying the damage linear density (expressed in vacancies ion −1 cm −1 ) by the implantation fluence (expressed in ions cm −2 ), thus yielding the volumetric vacancy density. This linear approximation is only suitable for low damage densities as it significantly over-estimates the real values at high damage densities where it does not account for nonlinear mechanisms such as self-annealing, ballistic annealing and defect interaction [44,51]. Accordingly, we have employed a more realistic description of the damage process as described below. Nonetheless, it is worth noting that the linear approximation gives a density of ∼2 × 10 23 vacancies cm −3 at the ion end of range, indicating that the sample has been implanted at a damage level that significantly exceeds the amorphization threshold in this region [44,[52][53][54][55][56]. The ion beam is characterized by a Gaussian point spread function. Therefore, to achieve a sharper edge between the irradiated and unirradiated areas, the implanted area was masked on its right side with the cleaved edge of a monocrystalline Si wafer. As shown in the inset of figure 2, the masked edge exhibits a sharper transition from the implanted to the unimplanted area although we might expect that there will also be scattering beneath the implant edge [57]. By taking into account the spatial resolution of the imaging technique, we estimated the sharpness of the edge as being better than ∼1 µm, which corresponds to the spatial resolution of the PL mapping techniques employed to characterize the sample (see below). As reported in figure 2, we will adopt a spatial reference convention so that the xand y-axes lie on the sample surface and correspond respectively to the directions perpendicular and parallel to the edge of the implanted region while the z-axis extends perpendicularly to the sample surface along the depth direction. The blue arrow in figure 2 highlights the direction of the PL scan to locate the 'strip' sector (reported below, see figure 8).
After ion implantation, the sample was thermally annealed with the purpose of forming a well-defined graphitic layer that would exert a stress field in the surrounding regions. The substrate was annealed for 1 h at a temperature of 800 • C, which is believed to be optimal for the formation of NV − aggregates [40], in a forming gas (4% H 2 in argon) ambient; particular care was taken to control the annealing atmosphere to avoid oxygen contamination which would have resulted in surface etching of the diamond structure.
PL mapping at micrometer resolution was performed on a Renishaw RM1000 micro-Raman spectrometer on regions surrounding the implanted area. The λ = 514.5 nm emission of an argon laser was employed as the excitation source after having been focused on the sample surface to a micrometer-sized spot by means of standard microscope optics with around 0.5-1 mW incident on the sample. The collection optics were configured for Raman backscattering measurements and comprised a notch filter for incident scattered light removal. The system was also equipped with computer-controlled precision xy stages to translate the sample under examination for two-dimensional mapping with ∼0.5 µm spatial resolution.
The spectrometer was equipped with a long-working-distance 50 × objective lens (0.7 NA) and a stage mounted cryostat allowing the acquisition micro-PL spectra and maps at liquid nitrogen temperature (T = 80 K). The use of a single grating with either 1200 or 1800 lines mm −1 coupled with a CCD array allowed the fast acquisition of maps with a spectral resolution that was adequate for PL measurements over a broad spectral range (500-800 nm). For both gratings, the spectral resolution of the system was evaluated as better than ∼0.17 THz (i.e. ∼0.2 nm) at a wavelength of 551 nm, by measuring the full-width at half-maximum of the firstorder Raman transition from the same diamond sample (i.e. ∼6 cm −1 ), thus fully adequate to measure the spectral features of interest (see below).
Using this configuration and the standard aperture (50 µm) of the confocal acquisition system, the probed depth inside the sample was estimated to be ∼1-2 µm. This implies that in the following analysis, the acquired signal should be considered as the average of the distribution of optical centers in the depth direction.
High-resolution micro-Raman mapping measurements were performed on a Dilor XY microspectrometer. The basic instrumental configurations (excitation wavelength, co-axial focusing/collection system, sample scanning system) are similar to that described for the PL mapping with the following notable differences: the sample objective lens was of a high magnification (80×), measurements were performed at room temperature and spectra were acquired in a triple-additive grating mode allowing for a higher spectral resolution over a narrower spectral range. The spectral width of the first-order Raman line in a high-quality diamond crystal was measured as 1.8 cm −1 . The laser power incident on the sample after focusing was ∼1.5 mW while confocal parameters determined a probe depth which was estimated to be within the ∼2-3 µm range.
Profilometry measurements were performed on the implanted sample both before and after thermal annealing to determine the swelling at each stage of sample processing. Profiles were obtained using an Ambios XP stylus profiler at a speed of 0.01 mm s −1 , force = 0.05 mg and filter setting 2. For each run the system was calibrated using a 186.2 nm vertical standard. A typical surface roughness of the sample before implantation was ∼5 nm.

Numerical simulations
To evaluate the stress dependence of the ZPL emission splitting and broadening of luminescent centers, numerical FEM simulations were performed to assess the stress fields established in diamond upon ion implantation. The stress derives from a decrease in mass-density and mechanical stiffness occurring in the material as a consequence of the ion damage to the crystal structure and resulting volume expansion of the amorphized regions. This volume variation is proportional to the mass-density variation and is partially inhibited by the mechanical reaction of the surrounding undamaged material. Due to the non-trivial geometry of the problem, FEM simulations are best suited to adequately model this constrained volume expansion.
The mass-density change in the amorphized region is not uniform, since it arises from the damage depth profile in all three dimensions (i.e. the depth profile in figure 3 convolved with the implantation aperture). To correctly estimate the mass-density (and hence volume) change, a recently introduced phenomenological model was adopted to describe the damage process in the crystal lattice at high implantation fluences [51]. The model is inspired by the work reported in [58] and takes into account the concentration of ion-induced vacancies with a simple linear approximation for the probability for a newly created vacancy to recombine with a self-interstitial: where P rec is the recombination probability at a given depth z and implantation fluence F and ρ V is the vacancy density in the material at depth z. α is an empirical parameter accounting for the defect recombination probability, i.e. it represents the saturation vacancy density in correspondence of which the recombination probability approaches 1 and therefore no further vacancies are induced in the structure. By solving the associated differential equation, we obtain an exponential saturation of the vacancy density for high implantation fluences. This relationship can be extended to the density variation in the damaged material as follows: where ρ d = 3.52 g cm −3 and ρ aC = 2.1 g cm −3 [44] are respectively the mass densities of pristine diamond and fully amorphized carbon and λ(z) is the linear density of induced vacancies per incoming ion as derived from SRIM. Starting from SRIM simulations for the 2 MeV He implantation, the depth profile of the mass-density variation in the as-implanted sample can therefore be derived (as shown in figure 3). Consistent with the approach adopted with the modeling of the mass-density variation, we assume that the mechanical properties of the asimplanted diamond (i.e. Young's modulus E(z) and Poisson's ratio ν(z)) also vary between the corresponding values of diamond (E d = 1220 GPa, ν = 0.20) and amorphous carbon (E aC = 21.38 GPa, ν = 0.45) with the same trend reported in (2) for the vacancy density. Following this approach, we obtained the three-dimensional distributions of mass-density and mechanical properties in the as-implanted sample, which we subsequently used as input in the FEM model to estimate the deformation and stress fields established in the sample [51].
To model the stress field after thermal annealing, for simplicity we invoke a simple binary approximation for the material composition. For regions where the vacancy density is below the graphitization threshold (D C ), we assume the lattice recovers to a crystalline structure. Conversely when the vacancy density exceeds D C , we assume that the material is fully converted to graphite. Hence, we modeled the depth profiles of the material properties (density, mechanical parameters) as step-like functions varying between values relevant to diamond and graphite (ρ g = 2.1 g cm −3 , E g = 10 GPa, ν g = 0.31), as shown in figure 3, where the dimensions of the graphitic region are set by the extent of the critically damaged region. With the annealed sample, the obtained three-dimensional mass-density distribution and spatially varying mechanical properties are input into the FEM numerical simulations, to estimate the deformation and stress fields established in the sample [59].
We note that our models include only two free empirical parameters, which are important at different stages of the processing: before annealing, α empirically accounts for the probability of an ion-induced vacancy/interstitial pair recombining with existing defects in the material, while for the post-annealing case D C is the critical density of vacancies above which the damaged structure converts to graphite upon thermal annealing. Note that while in principle these numbers should be known constants, in practice there is sufficient ambiguity in the literature as to their values, and hence treating them as parameters is appropriate. The other parameters, such as the mass-density and the mechanical properties of diamond, amorphized carbon and graphite, were taken from the accepted literature values [44,51,59].
The three-dimensional FEM model was adopted to calculate all stress components in the unimplanted regions close to the edge of the implanted area that result from the constrained expansion of the nearby buried damaged layer. A free mesh with tetrahedral elements was adopted and refined in the region of interest (i.e. the edge of the implanted region). In both cases (i.e. before and after thermal annealing) the non-uniform constrained expansion of the implanted crystal was directly determined by the local variations of mass-density and mechanical parameters and fitting applied to α and D C .

Results and discussion
We first report on the damage-induced swelling of diamond. As mentioned above, MeV ion damage induces a significant distortion of the crystal structure with accompanying swelling and induced pressure in the material [58,[60][61][62][63][64]. Figure 4 shows profilometry scans of the implanted sample before and after annealing, taken across the masked edge. Thermal annealing has the effect of increasing the swelling since the buried layer has been damaged beyond D C , and therefore the expansion of the buried layer dominates over the partial recovery (i.e. the compaction) of the upper cap layer. This result is consistent with previous reports on deep implantations [63]. Numerical FEM results, also shown in figure 4, are in good agreement with the experimental measurements for both the as-implanted and annealed samples, with best-fitting values of α = 4.4 × 10 22 vacancies cm −3 and D C = 2 × 10 22 vacancies cm −3 , respectively, in satisfactory agreement with the most recent reports [44,59]. The excellent agreement achieved between the experimental and numerical surface swelling estimations provides further confidence in the predictions of internal stress in the implanted diamond.
Shown in figure 5 are the simulation results for the post-annealing case, which is investigated in greater detail using the data from the PL mapping. In figure 5(a) we report a two-dimensional map of the four computed non-zero stress components (three principal and one shear) in a length/thickness (xz) cross-section of the specimen superimposed on the deformed shape of the specimen (where z displacements are multiplied by a factor of 2 to highlight the swelling effect). From figure 5(a), there is a strong compressive (i.e. positive, according to the sign convention adopted here) stress at the edge of the buried graphitic layer (highlighted in red) due to its constrained expansion. This stress produces a non-negligible effect on the investigated region of the sample near the edge of the implanted region (highlighted in blue), where both compressive and tensile stresses are observed at different depths and distances from the edge of the implanted area, due to the effect of the nearby swelling region. Since the Raman spectral features display a significant and well-established dependence on stress for diamond [65], the FEM results derived from the fitting of profilometry measurements were compared with those obtained from the micro-mapping of the first-order diamond Raman line. The shift and broadening of the peak was mapped along a linear scan extending perpendicularly from the edge of the masked side of the implanted area along the x direction (as shown schematically in the inset of figure 6(a)). From figure 6(a), the Raman peak measured in the stressed diamond region adjacent to the implanted region undergoes broadening and a shift determined by the varying stress fields. In particular, the evolution of the peak shift indicates that overall tensile stresses across the sample depth are established between distances 0 and ∼10 µm, while compressive components in the depth-averaged stresses dominate at distances >10 µm from the edge of the implanted region, until asymptotic zero-shift values are obtained at >30 µm. The dependence of the Raman peak shift ω on the hydrostatic stress σ h in diamond is given by the following empirical relation [65]: where a = 2.83 cm −1 GPa −1 and b = −3.65 × 10 -3 cm −1 GPa −2 . Through rearrangement, we derive the hydrostatic stress relation corresponding to the measured peak shift (as reported in figure 6(b)) together with the hydrostatic stress predicted through FEM simulations from the combination of the principal stress components (xx, yy and zz). In contrast to the PL measurements, the Raman data were estimated to originate from a confocal depth of ∼3 µm from the sample surface; therefore, FEM results were evaluated across this depth and then averaged. In the Raman measurements, the compressive components arising from deep regions are averaged with tensile components generated from shallow regions, thus resulting in a rather complex evolution of the measured 'mean stress'. This is reflected in the non-monotonic variation of the Raman-derived hydrostatic stress as a function of distance from the implanted region. As shown in figure 6(b), the agreement between the two datasets is satisfactory from an order-of-magnitude point of view, while the discrepancies between the experimental and numerical estimations of the hydrostatic stress are particularly pronounced where the shear xz stress is not negligible. This is understandable if we consider that (3) is strictly valid only when the principal stresses are applied to the crystal in an exclusively hydrostatic regime, while no empirical coefficients analogous to a and b are reported in the literature for the shear stress (xz). Nonetheless, the transition from overall (i.e. averaged between 0 and 3 µm) tensile stresses to overall compressive stresses at increasing distances from the edge of the implanted area can be intuitively explained if it is considered that at small distances stresses are mainly tensile at the surface and mainly compressive at depths of 2-3 µm (see figure 5(a)), while at larger distances stresses are more uniformly distributed through the depth and compressive effects dominate overall. The observed peak broadening can be ascribed to the non-uniform stress state, which removes the triple degeneracy of the first-order mode resulting in the splitting of the line into three components that could not be spectrally resolved. Scanning PL mapping of the ZPL NV − emission was performed on the as-implanted sample for the purpose of investigating the properties of the fraction of active luminescent centers that formed during the implantation process (i.e. before the thermal annealing). In particular, regions across the masked and unmasked edges of the implanted area were explored with linear scans in the main growth sector of the crystal as shown in figure 7. In all plots reported in figure 7, the black arrows indicate the scan direction consistent with the arrows in the inset pictures of figures 7(c)-(f). Therefore in all cases, the scans proceed away from the masked and unmasked edges of the implanted area. Figures 7(a) and (b) show PL spectra recorded along linear scans at 1 µm steps from the masked and unmasked edges, respectively. The spectra are normalized in intensity and displaced along the vertical axis to allow direct observation of the evolution of the peak shape. While the ZPL NV − emission shape remains substantially unchanged in the linear scan along the masked edge, a significant broadening is observed when scanning across the unmasked edge. The broadening increases as the probed point gets closer to the implanted region although a splitting of the λ = 637 nm emission into different spectral components is not clearly distinguishable (apart from a shoulder at λ = 635.5 nm). The different observed broadening behaviors in the masked and unmasked edges of the implanted area are attributed to the effect of stray implanted ions at the unmasked side of the implanted region. due to accidental surface tilting. Here we also observe significant differences between scans performed along the masked and unmasked edges. Firstly, the normalized PL intensity is much higher along the unmasked edge with respect to the masked edge while the two values tend to coincide asymptotically at larger distances from the implanted area. Once again we may attribute this observation to a larger concentration of stray ions implanted along the unmasked edge. Secondly, while the normalized PL intensity continues to increase with increasing distance from the implanted region across the masked edge, we observe a non-monotonic trend for the unmasked edge with a pronounced peak at 5-7 µm from the implanted region. This observation may be interpreted using the intensity of PL emission arising from defects in the diamond structure and also follow a non-monotonic trend. At low damage levels, the active PL centers are progressively created while at higher damage densities, the formation of optically active point defects is increasingly inhibited by the formation of more complex defects and ultimately by the formation of a continuous amorphous network [14]. This effect is visible at the unmasked edge where the spatial variation of the ion-induced defect density is broader. Along the masked edge, we only observe the native concentration of active NV − centers which are progressively quenched as they approach the implanted region. As shown in figures 7(e) and (f), the trend reported in figure 7(c) is monotonic on both the masked and unmasked edges as the scan progresses along the low-nitrogen sector highlighted in figure 2. Moreover, in the low-nitrogen 'strip' sector we measured similarly low PL intensities from the two edges of the implanted area.
Scanning PL mapping of the ZPL NV − emission was performed under the same experimental conditions on the annealed sample with the purpose of investigating the luminescent centers formed during the thermally activated process of recombination of nitrogen impurities and ion-induced vacancies in close spatial proximity to the heavily stressed implanted region. In this case, the linear scans were performed at 0.5 µm spatial steps. Figure 8 shows regions of the 'main' sector surrounding the implanted region, mapped using 60 µm long linear scans (#1-4). In figure 8, the edges of the implanted region are highlighted (red dashed line) while the coordinates of the starting points of the linear scans with respect to the respective corners of the implanted area are reported in blue. Figures 9(a)-(d) report the normalized PL spectra collected along the four scans reported in figure 8; in each plot, the corresponding scan direction is indicated as a black arrow. As observed in figure 7, the spectra are normalized in intensity and displaced along the vertical axis, to allow a direct observation of the evolution of the peak shape. As expected, after the thermal annealing of the sample, a significant increase of the NV − emission intensity was observed. We attribute this effect to PL-active centers created by straggling ions during the implantation process [57]. Moreover, the splitting of the NV − ZPL emission into up to four different spectral components is more clearly distinguishable with respect to the as-implanted sample. The position of the different components of the split NV − ZPL emission was evaluated by differentiating the spectra reported in figures 9(a)-(d) and identifying the position of the peak maxima in frequency units. The plots of figures 9(e)-(h) report the respective maxima as a function of position. The different spectral components induced by the strong stress fields located around the implanted region are distributed across a spectral range of ∼4 THz around the position of the non-perturbed NV − ZPL line (470.9 THz, corresponding to λ = 637 nm). The PL data also follow the local stress tensor and vary spatially due to the implant. The identification of the spectral components is clearer for the scans performed around the masked edge of the irradiated area (scans #1 and #2 reported respectively in figures 9(a) and (b)) where the transition between implanted and unimplanted regions is sharper. For scans #3 and #4 (figures 9(c) and (d)), peaks are broader and the identification of spectral components is less straightforward. For this reason, only results obtained from scans #1 and #2 will be commented on further.
Firstly, we observe how the non-perturbed line at 470.9 THz increasingly splits into three main components as the probed spot approaches the corner of the implanted region reaching a position located at 45 • with respect to the corner of the implanted area (i.e. the point of intersection of scans #1 and #2, at distance d = 20 µm with respect to their respective starting points). At this location, the shift in the three components with respect to the nonperturbed line is −0.5, +0.3 and +1.0 THz. As the scan progresses from the intersection point and approaches the lateral sides of the implanted area, a fourth component splits from the scans are performed along the unmasked edges of the implanted area, where we expect broader and more intense emission. A clearer identification of the spectral components of the ZPL NV − emission can be obtained by PL mapping the regions near the implanted area while moving along the 'strip' sector highlighted in figure 2. In this sector, the extremely low concentration of nitrogen impurities allows the measurement of much narrower spectral features, as shown in figure 10, where we report a linear scan across the 'strip' sector and far away from the implanted area. As with previous plots, the plot in figure 10(a) shows the spectra normalized in intensity and displaced along the vertical axis to allow a direct observation of the narrowing of the NV − emission when crossing the low-N sector. As reported in figure 10(b), the NV − emission intensity is also significantly reduced across the low-N sector. Figure 11 shows the result of PL linear mapping starting at the masked edge of the implanted area, and moving at 1 µm steps away from the implanted area while remaining within the low-nitrogen 'edge' sector (as shown schematically in the inset). As for previous images, the black arrows in the plots and in the insets indicate the scan direction. Figure 11(a) reports normalized PL spectra while figure 11(b) shows the variation of the spectral positions of the components of the NV − ZPL line in frequency units as a function of distance from the implanted area. As shown more clearly in a spectrum acquired for a longer integration time at a distance of ∼5 µm from the edge of the implanted area ( figure 12(a)), the narrower spectral components of the split NV − emission are clearly distinguishable and can be suitably fitted with a Lorentzian line function. Table 1 reports the outcomes of the fitting procedure. The splitting into up to five spectral components monotonically increases as the distance from the heavily stressed implanted area decreases, starting at about 15 µm from the edge and covering a range of ∼4.8 THz, i.e. a broader splitting range with respect to what was observed in the 'main' sector (see figure 9 for comparison). As a comparison, in figure 12(b) we report the PL spectrum acquired from a region at a distance of >20 µm, which we can assume to be unaffected by the stress field. As expected, a sharp single emission line is observed whose Lorentzian fitting parameters are also reported in table 1. With regard to the peaks linewidths, it is worth remarking that all PL spectra reported in this work represent the outcome of an ensemble measurement carried over a large number of active NV − centers included in the probing volume of the technique. In this sense, the significant broadening of the stress-split spectral components with respect to the zero-stress case is attributed to the averaging of PL signals from centers experiencing different local environment, i.e. orientation with respect to the lattice, local stress and surrounding defect/impurities.  [12]. In all plots, the black arrows indicate the scan direction, consistent with that reported in the inset of (b). and |EY states) and the applied stress [12]. If we apply such coefficients to the values of the different components of the complex stress field obtained from the numerical modeling (see section 3), it is possible to derive trends of the energy shifts for the different spectral components of NV − emission versus the distance from the implanted region. Table 2 reports the values of the coefficients, after conversion from energy to frequency units. The resulting trends, as reported in the continuous line plots in figure 11(b), are satisfactorily compatible with experimental data, particularly considering that the numerical model of the stress field has been derived under significant assumptions. These include (i) the same functional dependence of all mechanical properties with vacancy density, based on the outlined phenomenological model, and (ii) a simplified geometry of the buried graphite layer with respect to the real experimental situation. Moreover, the deconvolution of the multiple spectral components arising from stress fields along different 100 , 110 and 111 directions is not as unequivocal as in uniaxial stress experiments. Despite these assumptions, we observe that the stress 100 components in the x and y directions (i.e. where the stress field is more intense) can account for the observed trends. In the case of the |E X component arising from stress in the 100 direction (red curve in figure 11(b)), we have no experimental data at small distances from the implanted area where the numerical results would predict much larger shifts. Finally, we tentatively attribute the large shift at low frequency to a  shift of the |EY state induced by a stress field in the 101 direction. No unequivocal attribution could be found to relate the measured shifts to stresses in the 111 direction. Given the complexity of the observed effect, a fully detailed quantitative analysis of the observed trend is beyond the scope of this work and will be carried out in future systematic investigations.

Conclusions
We report the systematic investigation of the splitting of NV − photoluminescent centers induced in diamond by the stress field caused by MeV ion implantation. Engineering the strain field in implanted diamonds opens up new opportunities for NV − ensemble-based devices, by spectrally separating the inequivalent orientations and providing permanent shifts in the spectral features affording some degree of tuning. Ion-beam-induced stress fields are caused by the formation of the stable amorphous or graphitic phase in the highly damaged, sub-superficial region, respectively before and after thermal annealing. The ion-induced density variation creates complex internal stress fields in the unimplanted surrounding regions, resulting in the splitting of the NV − emission into different spectral components which have a range of ∼4.8 THz. The employment of implantation masks proved to be beneficial to improve the sharpness of the edges of the implanted areas. This allowed the definition of well-defined boundary geometries and in addition to the PL mapping of the NV − spectral features with micrometric spatial resolution in low-nitrogen sectors in the diamond crystal. Consistent with modeling in this work, the stress fields occurring in the implanted and annealed sample were simulated with FEM numerical codes by adopting a simple semianalytical model developed in previous work [51,58]. The free parameters in the model were determined by optimizing the consistency with surface swelling data while the model was further validated by testing its consistency with experimental Raman measurements.
The measured stress distributions at the edges of the implanted region were coupled with known coefficients linking static stress and NV − sub-level shifts [12], and the resulting values of the predicted line splitting proved to be in satisfactory agreement with the experimental data. A fully unequivocal attribution could not be established at this stage, and it will be the subject of future investigations.
Our results provide the first evidence for a new form of engineering of the properties of optical centers in diamond, namely implantation-induced strain engineering. In future, it should be possible to affect a particular strain field for a certain task by designing the appropriate implantation strategy for the optical response. This is significant as previous studies on applied strain in diamond have typically used external controls, such as diamond anvils [12], especially to reach the very high strains reported here. Such external control is impractical for most applications as the infrastructure required to maintain the strains is not compatible with the push for integrated devices.
The results obtained from the specific growth sectors of the sample show that the stressinduced splitting of the NV − emission can be measured with higher spectral resolution at lower nitrogen concentrations, thus indicating that future measurements of spin coherence, lifetime and polarization are possible at the single-center level in stress-engineered high-purity CVD samples, with significant potential applications in QIP.