Edge effects on optically detected magnetic resonance of vacancy defects in hexagonal boron nitride

The chemical and structural nature of defects responsible for quantum emission in hexagonal boron nitride (h-BN) remain unknown. Optically detected magnetic resonance (ODMR) measured from these defects was reported in two recent papers. In one case, the ODMR was tentatively attributed to the negatively charged boron vacancy, VB−\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_{\mathrm{B}}^ -$$\end{document}. Here we show how the key optical and magnetic properties vary with location within the bulk and along edges of h-BN sheets for this and the negatively charged nitrogen vacancy, VN−\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_{\mathrm{N}}^ -$$\end{document}. Sign changes of the axial zero-field interaction parameter D are predicted, as well interchange of singlet and triplet ground states. Based on the latest experimental information, we assign the observed ODMR signal to bulk VB−\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_{\mathrm{B}}^ -$$\end{document}. The other observed ODMR has some features reminiscent of our calculations for VN−\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_{\mathrm{N}}^ -$$\end{document} edge defects. Crystal defects in 2D materials play an important role in their physical properties and it important to understand the underlying physics to maximise certain characteristics. Here, the authors calculate how point defects in h-BN impact on its optical and magnetic properties and the variations which occur due to the location of the defect from bulk to edge.

T he discovery [1][2][3][4] in 2016 of single-photon emission (SPE) from defects in hexagonal boron nitride (h-BN) has inspired significant experimental and theoretical research in an effort to understand the chemical origin of this emission. Excitement stems from the possibility that such SPE could be harnessed by nanophotonic [5][6][7][8] industries to deliver applications such as quantum key distribution. Their exploitation, however, requires understanding of the chemical nature of the colour centres in h-BN and the ability to control and tune their structural arrangements and spectroscopic properties. Whilst various colour centres have been proposed and examined computationally as potential sources of SPEs in h-BN [9][10][11][12][13][14][15] , none have been firmly identified 13 . In parallel to this work, defect centres in h-BN have also been demonstrated to produce strong paramagnetic effects, and for these chemical assignments have indeed been developed 13 .
Recent advances, reported in two preliminary publications 16,17 , depict the observation of room-temperature optically detected magnetic resonance (ODMR) from h-BN, a phenomenon involving both measurement of ground-state magnetic properties and photoluminescence (PL) on the same defect. This not only makes it much easier to identify the chemical nature of a colour centre, but also suggests possible future applications in quantum information and other technologies.
In one preliminary ODMR report, by Gottscholl et al. 16 , based on the observed magnetic properties, the defect was tentatively assigned to the V À B defect, a negatively charged defect 9,11 , in which one boron atom is missing from the h-BN lattice. For this, a zero-field splitting (ZFS) parameter D of −3.5 GHz was reported for a triplet defect ground state, with a triplet ground state having previously been predicted for V À B 9,11 . Nevertheless, our initial calculations for a bulk V À B defect, reported herein, predicted a value of 3.32 GHz, in stark contrast to the reported value of −3.5 GHz. This inspired the remainder of the work presented herein, where we consider how the properties of defects may vary as the defect location is varied from bulk to some layer edge. This yielded another of our principle results, that the sign of D can change according to defect location. A relevant feature of the experiments is that the defects observed were located throughout the sample but were not native to the h-BN, instead produced following exposure to beams of neutrons, Li, or Ga atoms. Hence the bulk regions of the h-BN may sustain considerable damage, making the detailed atomic structure difficult to discern. This makes bulk versus edge differences all the more relevant; in native h-BN, defects located at both bulk 18,19 and edge 20,21 locations have been observed, as well as at structural defects in the basal plane of the h-BN 22 . Resolution of key issues appears imminent, however, as in the published version of the Gottscholl et al. 23 , the sign of D is changed to be positive. Another computational study confirms the large positive computed value of the zero-field splitting for V À B defects 16 . In the other preliminary ODMR report by Chejanovsky et al. 17 the ZFS parameter D was too small to measure (magnitude < 4 MHz) and, the hyperfine coupling constant (hfc) was determined to be around 10 MHz with an angular dependence consistent with an unpaired electron in a π orbital. In the report, the observed ODMR was tentatively associated with a substitutional carbon impurity or other low mass atom.
The most commonly applied computational approach to predicting the optical and magnetic properties of defects is densityfunctional theory (DFT) 13 . The prediction of optical spectroscopic properties, such as the zero phonon line (ZPL) energy, by DFT is very difficult as most defect states are open-shell in nature, displaying multi-reference characteristics, whereas DFT is intrinsically a single-reference approach 24 . Time-dependent DFT (TDDFT) can be reliable for defect excited states, provided that it is performed based on a well-represented reference state, as indeed is the case for the triplet manifold of V À B 13 . Nevertheless, DFT is typically quite accurate when it comes to the prediction of ground-state magnetic properties 13 , for example using the Heyd-Scuseria-Ernzerhof hybrid functional (HSE06) 25,26 . Indeed, DFT predictions of the hfc and ZFS properties of the nitrogen vacancy (NV) centre in diamond 27,28 were important for the identification of its chemical structure. Further, previously observed paramagnetic signals from h-BN defects have also been assigned by DFT as the neutral nitrogen (V N ) and boron (V B ) vacancies 29,30 . We therefore expect that appropriate combinations of DFT and TDDFT calculations will be reliable for the identification of defects associated with the observed ODMR.
Calculations pertaining to the Gottscholl et al. 16 ODMR observations have become available recently in the literature 31 . This also uses DFT to predict magnetic properties of bulk V À B defects, a critical novel conclusion being that strain-induced distortions away from full three-fold symmetry can interpret the observed off-axis ZFS parameter of E = 50 MHz. They evaluated excited-state energies using an ab initio multi-configurational approach, empirically modified to use DFT orbitals and solved using density-matrix renormalisation-group (DMRG) technologies.
The current works explores other possible explanations for the initially reported ZFS value reported by Gottscholl et al. 16 , as well as to develop understanding of bulk versus edge locality. We have also considered properties of the V À N defect which like V À B , has native 3-fold local point-group symmetry, D 3h . Our results suggest that the ODMR could instead be associated with edge-located V À N defects. As alternative defects, we also consider a defect involving an internal boron dangling bond 14 and the C þ N and C À N defects in which either a positively or negatively charged carbon replaces a nitrogen, respectively. We employ the CAM-B3LYP density functional [32][33][34] for this purpose, this calibrated, rangecorrected hybrid density functional being an example of an entrylevel method for the prediction of the spectroscopy of defects in molecular solids 24 . State energies reported by Ivády et al. 31 are of the same order of magnitude as the present results, but some differ in basic nature to key ones identified herein. Firstprinciples-based approaches analogous to theirs are available that may prove useful for modelling defect spectroscopy [35][36][37][38][39] , but are not applied herein.

Results
Zero-field splitting parameters. To find possible interpretations of reported ODMR, we investigated many possible edge locations for V À B , presenting also analogous results obtained for V À N and, in addition, some results for also a proposed 14 boron dangling-bond site as well as for C þ N and C À N . Initial investigations were performed examining calculated spin densities only, as these are straightforward to calculate reliably. Periodic sheets in either 1 or 2 dimensions are used in these investigations, the 1-dimensional structures presenting infinite straight h-BN edges. Results are given in detail in Supplementary Figs. 1-8, with Table 1 presenting a summary of some critical features of the calculated triplet ground states: defect natures, inner-ring geometry and location, and calculated value for the ZFS splitting parameter D.
A primary result is the demonstration that the ZFS parameter D can be very sensitive to the location of defects within the h-BN layer. For Hence all observable properties need to be individually considered. Concerning D, for V À N , the calculated values are also widely varying, ranging from −0.97 to 1.3 GHz. These results stress the need for detailed understanding of geometrical defect structure as a basis for the interpretation of ODMR, especially pertinent for experiments in which the defects of interest are produced following irradiation with particle beams 16,23 .
Concerning the recently observed ODMR, we select for further investigation from amongst the systems listed in Table 1: V À B in bulk, for which the calculated D agrees with the current experimental intepretations 23 , V À B at the postulated Hpassivated nitrogen edge, for which the calculated D value agrees with the original experimental interpretation 16 , and the analogous structures for V À N , both for reasons of exploiting the analogy to understand basic principles as well as for the possible interpretation of the ODMR experiments that reported |D| < 0.004 GHz as the calculated bulk value of −0.046 GHz may be in agreement with this to within the accuracy of the calculations 40  Origin of the edge effects for V À B and V À N . Figure 1a plots the calculated ZFS parameter D as a function of defect location for V À B . The one-dimensional strip geometric model used for this work is shown in Fig. 1b, with the calculated spin distributions of the triplet ground states superimposed. The corresponding data and geometric model for V À N are shown in Fig. 2a, b. In these figures, the defect is located adjacent to the edge, but the model supports moving the defect in towards the centre of the strip, emulating bulk defects. In both cases, the ZFS quickly approaches the bulk value as the defect moves away from the edge, but for V À N the change is small and regular, whereas for V À B an abrupt change reverses the sign of D upon moving just one site location away from the edge.
To understand the very different calculated qualitative properties, we note that lattice strains are already known to influence the spectral properties of some SPEs 41 . We compute D as a function of hydrostatic pressure, for both defects V À B and V À N in bulk, as shown in Fig. 3. This hydrostatic pressure produces isotropic compressive and tensile strains around the defect centres, thus preserving the D 3h point group symmetry, with D found to vary linearly with strain. In our calculations, the strain dependent ZFS yields a slope of 7.9 MHz/GPa (for positive strain) for V À B . This number, in principle, can be compared to the experimentally observed variation in ZFS with temperature, which has been attributed to thermal expansion of the lattice, yielding a slope of −0.4 MHz/K, equivalent to a pressure coefficient of 2.9 MHz/GPa 16 . These results do not entirely agree, unlike good agreement found between analogous calculations and observed properties for the NV centre in diamond 28 . If the defect were located at a step edge, then thermal expansion would lead to slippage in the direction orthogonal to the edge and hence the sensitivity of the ZFS would be reduced; alternatively, the location of the defect within the bulk sufficiently close to other defects caused perhaps by exposure to atom beams could allow for strain mitigation without affecting other defect properties.
One may distinguish two contributions to the variation of ZFS as a function of pressure: purely geometrical changes around the defect centre and the variation of the defect's spin density. The former may be described using the compressed-orbital model 28 , according to which the ZFS becomes scaled by a geometrical factor (d/d 0 ) determined by the atomic relaxation induced by pressure, in the proximity of the defect; here d and d 0 are the nearest-neighbour distances under pressure and at equilibrium, respectively. As shown in Fig. 3, the compressed-orbital model describes well the calculated pressure dependence of ZFS for both V À B and V À N , showing a negligible contribution arising from spin density changes. This is consistent with the results for NV centre in diamond, where the change in ZFS is attributed predominantly to geometrical changes 28,42 .
We now discuss the microscopic origin of the calculated change in D from the bulk values of 3.32 GHz in V À B and −46 MHz in V À N to the edge-values of −3.5 GHz and −30 MHz, respectively ( Table 1). As the edge is approached, the nearestneighbour bond lengths inside the defect change (see Supplementary Fig. 1) from three equivalent 2.59 Å separations to one at 2.59 Å and two at 2.81 Å for V À B ; for V À N , these distances are 2.04, and 2.30 Å or 2.04 Å, respectively, embodying the lowering of local point-group symmetry from D 3h to C 2v as the edge is approached. To demonstrate that geometrical changes to the inner-ring of defect atoms are responsible for the changes in ZFS parameters for V À B , we consider (Table 1 and Supplementary  Fig. 3) an artificially constructed structure in which the inner-ring atoms are frozen at their geometry as calculated at the boron edge, re-imbedding them inside the 2D bulk material. No significant change is calculated as a result of this re-imbedding, indicating that geometry alone dominates its large edge effect. For V À N for which the edge effect is small, analogous calculations indicate that the nature of the re-imbedding also becomes relevant (Table 1 and Supplementary Fig. 4).
The observed 16,23 non-zero value of E for V À B is non-zero may also be attributed to the edge effect. Our calculated value of E for the defect at the edge is 80 MHz, which on an absolute scale is in a reasonable agreement with the experimental value of 50 MHz 16,24 . Alternatively, residual strain influencing bulk defects could also account for the observed non-zero value of E 31 .
Photoluminescence and electronic states of V À B . For V À B , we focus on electronic states within the triplet manifold as previous calculations have suggested 9,11 that it has a triplet ground state, a property that is most likely required to explain the observed magnetic ODMR properties. All energetics calculations are performed using CAM-B3LYP on 3-ring model compounds described in detail in Supplementary Table 1 with TDDFT used to depict excited states. The triplet ground state in the bulk is predicted to be ð1Þ 3 A 0 2 , in accordance with recent calculations 31 and expectations, with the defect displaying D 3h local point group symmetry [9][10][11][12] . Such a ground state is consistent with the basic observations of Gottscholl et al. 16,23 . Note that we use standard group-theoretical axis conventions throughout 43 . When the defect is located at an H-terminated boron edge, the point-group symmetry is lowered to C 2v , with the ground state then described as (1) 3 B 2 . Table 2 lists critical calculated properties of the triplet manifold for defects located in the bulk and at the considered edge. These include the adiabatic transition energies ΔE 0 , likely to be a good approximation for observed ZPL energies ΔE 00 ; vertical excitation energies ΔE A v , and the absorption and emission reorganisation energies λ A and λ E , respectively, that depict spectral bandwidths 13 . State symmetry for bulk defects is expressed in terms of both D 3h to C 2v symmetry 43 , with spontaneous symmetry lowering mandated by Jahn-Teller distortion for doubly degenerate states; all edge defects intrinsically have C 2v symmetry; V À N in the centre (bulk) is also unstable owing to warping of the plane to C 2 symmetry. The CAM-B3LYP TDDFT calculations are expected to be very reliable for the situation presented by V À B but typically underestimate isolated-layer transition energies by 0.3 ± 0.2 eV in situations relevant to V À N at the edge, increasing to 0.4 ± 0.3 eV for V À N at the centre 13 . No corrections are made for the embedding of the defect layer inside bulk h-BN, but cluster calculations show small change as the ring size is increased from 1 to 2 and ca. 0.03 eV in expansion from 2 to 3. The (1) 3 B 1 , (1) 3 A 2 states of the V À B defect located in the centre and edge, respectively, presumably form a transition-state on the (possibly distorted) Jahn-Teller potential-energy surface. In addition the latter state along with the (1) 3 A 2 state at the centre and (2)  states for the edge location are Franck-Condon allowed. The adsorption vertical oscillator strength for the (1) 1 A 2 state of V À N located at the boron zigzag edge is 0.028. When λ A þ λ E > E A v , an energetically accessible conical intersection with the ground state is present, making an ultrafast relaxation pathway for initial excitation back to the ground state as well as for relaxation to a potentially stable defect isomer. This is the case for the (1) 3 A 2 and (1) 3 B 1 states of the V À N defect located at the centre. The lowest vertical excitations from these two triplet states are of B 1 type, at 2.19 eV for (1) 3 A 2 and 1.91 eV for (1) 3 B 1 .
Details of the dominant orbital excitations contributing to the example system of bulk V À B at vertical excitation are presented in SI Section S3; this includes also a description of spincontamination effects, and embodies small contributions from h-BN valence-band excitations to the excited states, listing all five excited-state components of up to 3.5 eV in energy.
In the centre of an h-BN layer, its first triplet excited state is predicted to be (1) 3 E″, a state that undergoes a Jahn-Teller distortion to produce the (1) 3 A 2 state in C 2v symmetry at an adiabatic transition energy of 1.71 eV. This quantity is a good approximation for the ZPL energy and is in excellent agreement a b Fig. 3 Pressure dependence of the zero-field splitting parameter for negatively charged boron and nitrogen vacancies V À B and V À N . a Calculated pressure dependence of the axial component of the zero-field splitting parameter (ZFS), D, as a function of strain applied to the lattice for V À B using two different models. The green line shows the results of a complete self-consistent solution (with geometry relaxation) of the strained supercells of h-BN embedding a single defect. The purple line represents the ZFS calculated according to analytical compressed orbital model described in ref. 28 . In this model, only the distance of the spins is affected, while the change of the spin density distribution is neglected. b The same calculation for the V À N defect. Table 2 Time dependent density functional theory (TD-DFT) calculations using the CAM-B3LYP exchange-correlation functional for the negatively charged boron vacancy V À B and negatively charged nitrogen vacancy V À N .

Defect
In centre At boron zigzag edge A 3-ring model compound is used with the defect located at the centre or at the boron zigzag edge of the model. State symmetry for bulk defects are expressed in terms of both D3h to C2v symmetry. Calculated state properties are the adiabatic transition energy ΔE0 (in eV), the vertical excitation energy ΔE A v (in eV), the absorption and emission reorganisation energies λ A and λ E , respectively (in eV).
with observed 16 photoluminescence properties. The transition is formally forbidden in D 3h symmetry, which may explain the very weak observed emission as the implied long radiative lifetimes could permit other processes to operate. The transition becomes allowed out-of-plane polarised once the Jahn-Teller distortion occurs. Also, (1) 3 E″ appears to be very close in energy to ð1Þ 3 A 00 1 , a state that inherently supports allowed emission to the ground state (calculated transition moment 0.18 D). Of note too, ð1Þ 3 A 00 1 is predicted to undergo an avoided crossing with (1) 3 A 2 as the symmetry is lowered to C 2v .
A critical feature is that the calculated reorganisation energy in emission λ E of 0.29 eV is at least double the largest values that could be consistent with the observed spectral bandwidth 16,23 if simple approaches that ignore the Jahn-Teller effect are used to model the spectra. The observed spectra are noted as arriving from many emitters rather than single-photon sources, and it is possible that the observed bandwidths reflect inhomogeneity rather than intrinsic emitter spectral width. Improved experimental data is required so that the extent of the discrepancy between observation and calculation can be assessed. Possibilities that need to be investigated in future work include that the observed emission is significantly narrowed from naïve expectations owing to either the Jahn-Teller effect or the impact of acoustic phonons, and that the emission could arise from ð1Þ 3 A 00 1 , or else from within the singlet manifold 31 following intersystem crossing. Table 2 also shows changes to the calculated state properties for V À B located at a H-terminated boron edge. Mostly, the changes are small, yet the impact could be profound as, e.g., the (1) 3 A 2 and (1) 3 B 1 structures reverse in relative energy, dramatically changing the nature of the (now pseudo) Jahn-Teller distortion. Most significantly, as expected, the calculated emission reorganisation energies λ E increase significantly, becoming 0.53 eV for the lowest-energy triplet photoemission, depicting spectra very much broader than observed. Hence it is unlikely that the observed ODMR involves the triplet manifold of V À B located at this edge.
Photoluminescence and the electronic states of V À N . The V À N defect in the bulk has been predicted to have a closed-shell singlet ð1Þ 1 A 0 1 ground state [10][11][12]30 , and our calculations support this prediction. ODMR can be observed for systems with singlet ground states that allow for intersystem crossings (ISC) from initially excited singlet states to triplet states 44 , but this requires also significant buildup of population on the triplet state and hence unusual circumstances. The calculated photoluminescence properties of the singlet ground state are listed in Table 2 and indicate that it is highly unlikely that this could lead to any of the observed ODMR signals: these signals occur following excitation by light with under 2 eV energy per photon, whereas the lowestenergy singlet excited state is predicted to be located vertically at 3 eV. Further, the lowest excited state is predicted to undergo a huge Jahn-Teller distortion with a reorganisation energy of 1.19 eV, two orders of magnitude larger than observed spectral widths. It also manifests out-of-plane buckling and a conical intersection with the ground state that would facilitate extremely fast non-radiative energy loss.
Of note, however, is that our calculations predict that the lowest-energy triplet states are just 0.3-0.5 eV higher in energy than ð1Þ 1 A 0 1 , numbers of the order of likely errors 13 in the computational approach and hence the ground state may indeed be triplet in character. The lowest-energy transitions within the triplet manifold are predicted to be out-of-plane polarised at 2.19 eV for (1) 3 A 2 (Jahn-Teller transition states) and 1.91 eV for (1) 3 B 1 (Jahn-Teller local minima), consistent with observed PL properties. Table 2 also shows that moving the defect to an edge has a profound effect on the properties of the ground and excited states of V À N , making a firm prediction of a (1) 3 B 1 ground state. For the edge-located defect used in the calculations (see Table 2), the excitation energies within the triplet manifold now appear too large to account for the ODMR results, as also do the reorganisation energies. Nevertheless, the calculations raise the possibility that some other edge-related V À N defect could have the observed PL properties as they are established as being strongly site dependent. They show that relocating defects from the bulk to edges can have profound effects on electronic properties, analogous to the profound effects on magnetic properties predicted for V À B .

Discussion
Our calculations show that the optical and magnetic properties of defects may change considerably as a function of the location of the defect in regard to step edges in h-BN. For V À B , the calculated ZFS parameter D changed sign on going from edge to centre, with only small changes to the PL, whereas for V À N , D only changed slightly for the triplet state but the nature of the ground state changed from triplet to singlet. Defect location and defect properties are strongly correlated.
The calculations presented here for ground state spin properties and the PL spectra provide strong evidence that the recently observed ODMR/EPR spectra for h-BN by Gottscholl et al. 16,23 arises from the V À B defect. Our calculated values for the magnitude of the ZFS parameter D and hfc tensor A for bulk-located defects are in excellent agreement with those finally reported by Gottscholl et al. 23 . Nevertheless, ignoring possible spectral changes owing to the Jahn-Teller effect, the calculated reorganisation energies for emission within the triplet manifold appear to be too large to explain the narrow observed photoemission linewidths obtained for an ensemble of emitters, and further experiments and calculations are required to understand native spectral line widths, the significance of the Jahn-Teller effect, the role of acoustic phonons, and alternate emission possibilities.
Also, the calculations suggest that the ODMR observations of Chejanovsky et al. 17 could possibly be assigned to V À N defects located near step edges. At the edge, the calculations strongly suggest a triplet ground state, with a singlet ground state slightly preferred in the centre with PL properties similar to those observed. The ZFS parameter D varies between −46 and −30 MHz for location in the bulk and at the edge respectively and hfc tensor of A = [10.0, 10.0, 29.5] MHz. These results are in a reasonable agreement with the experimental observations, whereas one considered simple defect involving carbon, as was originally suggested as a possible centre, was found to be inconsistent. The relationship of these results to more widely obtained results for other SPE's 13 needs to be established.

Methods
The computational parameters and theoretical details are presented in detail in Supplementary Methods 25,26,[45][46][47][48] . Briefly, periodic DFT calculations are performed by the Vienna Ab Initio Simulation Package (VASP) 45,46 using the nonlocal Heyd-Scuseria-Ernzerhof hybrid functional (HSE06) 25,26 in replicated images representing h-BN monolayers, strips, and multi-layer slabs. For the strips, we assume that all h-BN edges are hydrogen terminated, with no access of hydrogens to partially fill defect vacancy sites, even if the defect is located adjacent to an edge. Magnetic properties are evaluated using standard methods 28,30,49,50 . Excited-state properties are evaluated on a 3-ring model compound (see Supplementary Data) representing the defects, evaluating transition energies, reorganisation energies, and oscillator strengths by TDDFT using Gaussian-16 51 . The CAM-B3LYP density functional 32-34 is used with the 6-31G* basis set 52 , this being a calibrated entry-level method for the evaluation of defect spectroscopic properties 24 .