A planar defect spin sensor in a two-dimensional material susceptible to strain and electric fields

The boron-vacancy spin defect ($\text{V}_\text{B}^{-}$) in hexagonal boron nitride (hBN) has a great potential as a quantum sensor in a two-dimensional material that can directly probe various external perturbations in atomic-scale proximity to the quantum sensing layer. Here, we apply first principles calculations to determine the coupling of the $\text{V}_\text{B}^{-}$ electronic spin to strain and electric fields. Our work unravels the interplay between local piezoelectric and elastic effects contributing to the final response to the electric fields. The theoretical predictions are then used to analyse optically detected magnetic resonance (ODMR) spectra recorded on hBN crystals containing different densities of $\text{V}_\text{B}^{-}$ centres. We prove that the orthorhombic zero-field splitting parameter results from local electric fields produced by surrounding charge defects. By providing calculations of the spin-strain and spin-electric field couplings, this work paves the way towards applications of $\text{V}_\text{B}^{-}$ centres for quantitative electric field imaging and quantum sensing under pressure.


I. INTRODUCTION
Electric and strain field effects are highly important ingredients in the tight control of solid-state spin defects for quantum technology applications.The accurate knowledge of the couplings to these perturbations can be harnessed either actively, e.g. for tuning the magneto-optical properties of the defect, or passively, e.g. for sensing its close environment.Particular attention is currently given to spin defects embedded in two-dimensional (2D) van der Waals materials [1][2][3][4].The main advantage of a 2D host is the proximity of the quantum defect to the surface, which improves both its optical efficiency [5] and sensing capabilities [4].One of these 2D materials is hexagonal boron nitride (hBN), which can be exfoliated down to the mono-layer limit without compromising its chemical stability [6,7].Owing to its large optical band gap (∼ 6.0 eV [8]), hBN hosts optically-active point defects in a wide range of wavelengths.Once isolated at the individual scale, defects in hBN were first studied as single photon emitters, providing narrow and tunable emission lines, high brightness, and perfect photostability [9][10][11][12][13][14][15][16][17].In addition, it was recently shown that the electron spin resonance frequencies of some defects can be inferred through optically detected magnetic resonance (ODMR) methods, providing a central resource for quantum sensing applications [18][19][20][21][22][23][24].
Under optical illumination, the V − B centre exhibits a very broad emission in the near-infrared [24].Although an experimental work has reported a zero-phonon-line (ZPL) at 773 nm [41], a recent theoretical study has claimed that the full emission is phonon-assisted since optical transitions between electronic states are forbidden [42].As a consequence, the optical emission is intrinsically dim and only dense ensembles of V − B centres could be detected in experiments so far.Theory further identified that the V − B centre has a spin triplet ground state (S = 1) with D 3h symmetry [26].In the absence of external perturbations, the ground state spin Hamiltonian only involves spin-spin interaction and reads as The axial zero-field-splitting (ZFS) parameter D separates the m S = 0 and m S = ±1 spin sublevels, the notation m S referring to the spin projection along the caxis (z) of the hBN crystal.Interestingly, ODMR spectra recorded on ensembles of V − B centres usually indicate two magnetic resonances at frequencies D ± E, with D ∼ 3.47 GHz and E ∼ 50 MHz [24].These measurements thus reveal an additional orthorhombic splitting (E) of the m S = ±1 levels, which is in contradiction with the D 3h symmetry of the V − B centre.While this splitting was originally attributed to strain effects in the hBN crystal leading to a reduced symmetry of the defect [24,31], a recent study suggested that it results instead from the interaction with a local electric field [43].However, no microscopic picture has been applied to interpret these results.More generally, there is an urgent need to understand the coupling of strain and electric fields to V − defect spins, which is an exemplary quantum sensor in a van der Waals material.
In this paper, we employ first principles simulations to determine the effect of external electric and strain fields on the ground state electronic spin structure of V − B centres in hBN.We show that in most of the experimental conditions the defect experiences fluctuating electric fields that leads to an orthorhombic splitting in ODMR spectra recorded at zero external magnetic field.We identify the microscopic origin of this effect and determine the spin-strain and spin-electric field coupling parameters.We then compare these results to experimental ODMR spectra recorded from neutron-irradiated hBN crystals containing different densities of V − B centres.Our study proves that the orthorhombic ZFS is caused by fluctuating electric fields and suggest that the coupling strengths of the V − B centres to strain and electric fields are comparable to those of the nitrogen-vacancy (NV − ) centre in diamond [44,45].

A. Defect structure and spin-spin interaction
The V − B centre in hBN is a single vacancy at the boron site in the negatively charged state possessing D 3h point symmetry (Fig. 1a).In the following, the spin properties of the V − B centre are calculated using a large hydrogenterminated mono-layer flake model (see Methods).The coordinate frame is defined as {x, y} span the plane of the hBN layer, with x and y along the zigzag and armchair directions of the hexagonal lattice, respectively.The Kohn-Sham levels in the spin-polarised calculations are plotted in Fig. 1b, showing nearby double degenerate e and e levels and non-degenerate a 1 and a 2 levels leaving two holes in the e orbital of the spin minority spin channel.This constitutes the 3 A 2 spin triplet ground state of the V − B centre, in agreement with previous calculations [25,26,46].Importantly, the spin density is dominantly localised in-plane on the three dangling bonds of the neighbouring nitrogen atoms (Fig. 1a), as recently confirmed experimentally by pulsed electron-nuclear double resonance techniques [47].
Starting with an unperturbed D 3h symmetry, the mono-layer flake model leads to an axial ZFS parameter D = 3263 MHz [see Eq. ( 1)], a value in fair agreement with the experimental results D ≈ 3470 MHz.In the Supplementary Note 1, we show that the mono-layer flake model well represents the bulk environment around the defect because the spin density is fully localised into the hBN plane.Although there is about 6% discrepancy in the absolute value of the D parameters obtained from our theoretical approach and the experimental data, the accuracy in the variation of the zero-field-splitting parameters upon external perturbations can be reliably obtained thanks to the cancellation of numerical errors in partial derivatives.In the next sections, we compute the coupling coefficients of the V − B centre to strain and electric fields.We assume that the strengths of theses fields are much smaller than the Coulomb interaction between the electrons, so that strain and electric fields can be considered as perturbations.We analyse how these interactions can lead to an orthorhombic ZFS parameter, that splits the m S = ±1 levels upon axial symmetry breaking.

B. Strain field interaction
We start the analysis of external perturbations with the strain field case.We construct the interaction Hamiltonian analogously to previous calculations performed on the NV − centre in diamond featuring C 3v symmetry [48].In the linear approximation, the bi-linear forms of the second rank tensor of the external strain are connected to the bi-linear forms of the Ŝ • Ŝ second rank tensor by scalar coupling coefficients g.We write the interaction in a symmetry-adapted form to retain only the linearly independent couplings, i.e. we construct the irreducible representation basis in D 3h symmetry from the components of the strain and spin tensors.These transform as quadratic spatial coordinates, e.g.zz and xx + yy both belong to the totally symmetric A 1 irreducible representation.The interaction Hamiltonian is totally symmetric, thus only the products of the same irreducible representations of the strain and spin are valid combinations.
According to the theory of invariants, the Hamiltonian describing the interaction of the V − B centre with strain can be formulated as where g i are the coupling coefficients between ε bi-linear forms of the strain matrix elements and Ŝ bi-linear forms of the spin projection operators reflecting a symmetryadapted form of A 1 × A 1 + E × E .Note that the out-of-plane strain components {ε zx , ε zy , ε zz } are not included in the Hamiltonian because we consider a singlelayer flake model.However, these couplings are expected to be minor for a bulk model as well since out-of-plane (z) distortions have a small impact on the in-plane hBN structure that contains the spin density of the V − B centre (see Supplementary Note 2).
To calculate the coupling coefficients g i , we specify uniform strain of a single element in the strain tensor in each calculation and apply it to our model (see Methods).We do so for each symmetrically non-equivalent components in Eq. ( 2), i.e. for ε xx and ε xy .We then calculate the resulting axial and orthorhombic ZFS parameters (see red data in Fig. 2a,b) and obtain the coupling coefficients g i as the partial derivatives where D ij are components of the ZFS D-matrix defined by alternate form of the interaction Ĥ = ŜD Ŝ.The negative sign of g 1 implies a larger (smaller) axial ZFS parameter D = 3D zz /2 for compressive (tensile) external strain, in line with the decrease (increase) of the distance between the localised spin density lobes.We now analyse how the microscopic structure of the defect affects the response to the strain.To this end, we carried out calculations for which we strained the lattice of the defective model without allowing any local ionic relaxation around the defect and then calculated the change in the ZFS parameters (blue points in Fig. 2a,b).After allowing local relaxation under the externally applied strain (red points in Fig. 2a,b), we can define the local geometry changes in the close vicinity of the defect.The local strain is defined from the deformation of the triangle spanned by the three neighbouring nitrogen atom positions.We plot this local strain as a function of the external ε xx strain in Fig. 2c.We identify a large b Orthorhombic ZFS parameter E as a function of the external strain.The cases of fixed and relaxed atomic positions are compared (blue and red points), as explained in the main text.c Local strain response for the externally applied strain.Blue dots show the one-to-one correspondence when the atomic positions are not allowed to relax.Red and green points show the enhanced local εxx strain and the activated local εyy strain around the defect, respectively, when its close vicinity is allowed to relax under the strain applied to the crystal.
increase in the same component (ε xx ) of the local strain after relaxation.Moreover, an additional ε yy component is also activated.We attribute this effect to the smaller local stiffness at the vacancy site originating from the broken bonds.During the relaxation of the defect under external strain, it is energetically favourable for the vacancy site to accommodate an increased strain compared to the externally applied one, thus lowering the strain enthalpy of the host material in its vicinity.Generally, the same strain enhancement is expected in defects with dangling or weaker bonds than the ones in its host material, making this type of qubits more sensitive strain sensors.The increased local strain after relaxation reflects in the increase of the ZFS strain coupling strength as well.Note that relaxation under external sheer strain (ε xy ) shows similar trend in the ZFS coupling strength.
Finally, we convert the strain-spin coupling strengths to stress couplings, as stress is usually more accessible and controllable in experiments than strain.Based on the experimentally available elastic parameters of C 11 = 811 GPa and C 12 = 168 GPa in bulk hBN [49], we calculate the h i stress coupling coefficients as where C ij elements are the second order elastic constants (stiffness tensor) in the Voigt notation.The relation between the strain (ε) and stress (σ) can be then formulated in linear elasticity as With the above conversion, the Hamiltonian takes a similar form to Eq. ( 2) using the stress coupling coefficients Importantly, our calculations already indicate that |h 1 | > |h 2 |, such that spin-stress coupling mostly leads to variations of the axial ZFS parameter D without inducing a significant orthorhombic splitting.

C. Electric field interaction
The Hamiltonian of the electric field coupling is constructed analogously to the strain coupling as where d ⊥ is the coupling coefficient between E linear forms of the electric field vector and Ŝ bi-linear forms of the spin projection operators reflecting a symmetryadapted form of E × E .We note that the parallel component of the electric field (E z ) with respect to the hBN plane normal vector belongs to the A 2 irreducible representation, and it cannot couple to either A 1 or E representations of the spin operators.As a result, the interaction with an electric field only introduces an orthorhombic splitting without inducing variations of the axial ZFS parameter D.
We apply homogeneous electric field to the defective hBN model.The coupling coefficient can be calculated as the partial derivative of the ZFS parameter (see Fig. 3a) We note that a recent experimental work has estimated the coupling strength of the E parameter to the electric field as d ⊥ = 35 Hz cm/V [43], which is in the same order of magnitude as our first-principle results.
We now analyse the geometry relaxation under the applied external electric field similarly as for the case of the strain perturbation discussed above.We can identify an accompanying piezo effect resulting in (ε xx − ε yy ) and (ε xy + ε yx ) additional strains for the applied E y and E x electric fields, respectively.The former calculation is shown in Fig. 3b.This effect enhances the coupling strength to the electric field by a factor of 1.5 (compare blue and red points in Fig. 3a).
The microscopic origin of the piezo enhancement effect can be explained by changes in electrostatic interactions in the vicinity of the defect under the applied external electric field.We obtain a significant change in the partial charge distribution of the strongly polarised bonds around the vacancy defect depending on the relative direction of the external electric field.This implies a redistribution of the electron density by the electric fields, which generates quantum mechanical forces on the ions.The origin of this effect can be detailed in the example of the vacancy-neighbour nitrogen atom in the y direction and its boron neighbours.Without the external electric field applied, the B-N bonds are polarised with an effective dipole moment pointing from N to B. The applied field in the positive (resp.negative) direction increases (resp.decreases) the original charge separation and the effective dipole moment.This results in an additional effective attractive (resp.repulsive) force along the B-N bond and an outward (resp.inward) relaxation of the nitrogen atom relative to the vacancy site.We identify this additional relaxation as the ε yy piezo response.Thus, the relaxation of the nitrogen atom is in the same direction as the electric field.This finally enhances the symmetrybreaking effect of the electric field and consequently its coupling to the E parameter.

D. ODMR simulation
In the experimental ODMR signal of the V − B centre, a small orthorhombic E splitting was reported despite its contradiction with the defect symmetry.Recently, the effect was attributed to the local electric field originating from nearby charged defects without providing arguments from microscopic picture of the defect [43].Firstly, we simulate the ZFS parameters of a single V − B defects affected by the presence of uniformly distributed random local strain or electric fields of given magnitudes.We sample 1000 random configurations for each given magnitude and visualise the corresponding m S = ±1 eigenvalues of the ZFS Hamiltonian as a scatter plot as a function of the perturbation magnitude.From the simulations shown in Fig. 4, we conclude that the E splitting can be exclusively attributed to the local electric field, while the effect of strain leads to a broadening of the signal.The sharp splitting originates from the large difference in the magnitudes of the couplings of the D and E parameters to the external fields.We can describe this with the ratio g D /g E of the coupling strengths.In the case of the electric field, the D parameter coupling is symmetry forbidden, g D /g E ≈ 0, leading to a sharp splitting of the magnetic resonances.In the case of the strain field, we calculate g D /g E = 7.38, which leads to a dominating shift in the ODMR signal with broadening of the integrated spectrum.We note that a similar effect was already observed and simulated for the NV − centre in diamond [50].
We now use our theoretical predictions to analyse experimental ODMR spectra recorded at zero external magnetic field.To this end, we rely on three hBN crystals (S1, S2 and S3) with different densities of V − B centres created by neutron irradiation.Their optical properties are analysed with a confocal microscope operating at room temperature under green laser illumination (see Methods).All crystals exhibit the characteristic emission spectrum of V − B centres with a broad spectral line in the near infrared.The relative density of V − B centres is estimated by recording the evolution of the PL signal with the optical excitation power [Fig.5a].In the considered power range, the PL signal increases linearly with a slope proportional to the density ρ of V − B centres.By analysing these data, we obtain ρ S2 ∼ 4.8 × ρ S1 and ρ S3 ∼ 8.5 × ρ S1 .
The ODMR spectra recorded at zero external magnetic field are shown in Fig. 5b.For the three crystals, we detect the two characteristic magnetic resonances of the V − B centre with frequencies at D ±E.Importantly, the D parameter is identical for all spectra while the E-splitting increases with the density of V − B centre.To understand these results, we rely on a microscopic charge model orig- inally introduced for NV − defects in diamond [50].We consider that negatively-charged V − B centres are associated to positively-charged defects in order to ensure charge neutrality of the hBN crystal.These charges produce a local electric field that depends on the charge density ρ c .The Hamiltonian describing the interaction of a central V − B centre with this electric field is constructed from Eqs. ( 1) and ( 13) using the calculated coupling coefficient d ⊥ .Moreover, we include the hyperfine splittings from the first neighbour 14 N nuclear spins with a coupling strength of 47 MHz [26,51].We keep this value fixed in the simulation according to the obtained negligible effect of the external electric field on the hyperfine parameters.We note that the hyperfine constants are, however, sensitive to the strain (see Supplementary Note 3).For the microscopic model of the charge environment, a specific number of elementary point charges are placed at the atomic sites of hBN in a simulation radius of 10 nm, corresponding to an average charge density ρ c in the simulation sphere.Their position vectors relative to the V − B defect at the origin of the sphere is sampled randomly from the standard normal distribution.The resulting total electric field at the origin is supplied to the x o 1 x a 9 x 9 p B o D h W Y B 3 5 Z x / w 7 A 2 J n 9 < / l a t e x i t > ⇢ c = 0.018(7) nm 3 < l a t e x i t s h a 1 _ b a s e 6 4 = " b k W g e y + S F x J q r r /  Hamiltonian as an effective external electric field.The obtained spectrum is collected for 10 4 different random charge configurations.We then optimise the charge density and the ODMR contrast to fit to the experimental spectra.
As shown in Fig. 5b, all ODMR spectra are well fitted by this microscopic charge model, leading to charge densities ρ c,S1 = 0.018(3) nm −3 , ρ c,S2 = 0.046(4) nm −3 and ρ c,S3 = 0.141(8) nm −3 .These results indicate that the density of charges increases linearly with the density of V − B centres in the hBN crystal (Fig. 5c), as expected from the microscopic charge model.

III. DISCUSSION
We calculated the coupling strengths of the ZFS parameters to the symmetry adapted components of strain and electric fields for the V − B defect in 2D hBN from first principles.In order to place these values in a broader context, we briefly compare them to those of the popular NV − defect in diamond.The calculated and exper-imental values of the stress coupling on the axial ZFS parameter of the NV − centre are −5.17MHz/GPa [48] and −4.4 MHz/GPa [44], respectively.These values are weaker than that obtained for the V − B centre h 1 = −19.6MHz/GPa.On the other hand, the experimental perpendicular electric field coupling d ⊥ was reported to be 17 Hz cm/V for NV − centres [45], which is close to the result of our calculation of about 21 Hz cm/V for the V − B centre in hBN.
For the strain perturbation, the coupling to the axial ZFS parameter D dominates, while it is forbidden by symmetry for electric field perturbations.Consequently, we identify the electric field perturbation as the source of the symmetry-breaking orthorhombic E splitting in the ODMR experiments.Our numerical simulations based on the calculated couplings can accurately model the experimental ODMR signal in the presence of electric field perturbations originating from point charges in the lattice sites of the bulk hBN host.Furthermore, we identify a correlation in the defect density created by neutron irradiation with the simulated density of point charges and ultimately with the effective E splitting in the ODMR signal.
Moreover, our DFT calculations reveal a local enhancement effect in the coupling strengths for both strain and electric fields.The former is specific to the vacancy nature of the defect.The latter is attributed to the local piezo effect, i.e. an accompanying strain response under the external electric field effect.In other words, the local piezo effect strengthens the response of the ZFS to the external electric field which is a significant effect of about 50%.We note that this effect goes to the opposite of one 3D qubit system, the divacancy qubit in 4H silicon carbide [52].In that case, the local piezo effect rather generates such quantum mechanical forces on the ions that compensate the effect of the external electric field and reduce the spin-electric field coupling by about 10%.
In summary, the spin-strain and spin-electric field couplings of V − B centre in hBN were determined.We proved that fluctuating electric fields around the V − B spins are at the origin of the orthorhombic splitting commonly observed in ODMR spectra of ensembles V − B centres at zero external magnetic field.We showed that the piezo effect is strong in the spin-electric field coupling parameters, which may be generalised for other vacancy-type defects in hBN.This work might guide future applications of V − B centres in hBN, such as quantum sensing under high pressure and electrometry.

A. Density functional theory calculations
The nano-flake model used in the calculation is a hydrogen terminated hBN mono-layer of up-to the fourth hexagon neighbours around the boron vacancy defect, consisting of 147 atoms in total.Its electronic structure was calculated by density functional theory (DFT) using the PBE0 functional [53] and split valence polarisation Karlsruhe basis set (def2-SVP) [54], as implemented in the ORCA code [55].The spin-spin ZFS tensor is calculated using the restricted spin density obtained from the singly occupied unrestricted natural orbitals [56].During the relaxation of the defect structure under external strain applied, the outshell hBN hexagons were fixed in the atomic positions that are associated the strain environment acting to the inner part of the model.

B. Numerical simulations
ZFS simulation of the defect under external strain and electric field perturbations is calculated in a home-built Python code.We iterate over the magnitude of the perturbation in 1000 steps for the simulation range.In each iteration, we adaptively sample (scaled by the magnitude) the random configurations of effective perturbation vectors ε = (ε xx , ε yy ) and E = (E x , E y ) for the normal strain and electric fields, respectively.The component are sampled from uniform random distribution and the vector is normalised.In the ODMR simulation, we sample a uniform random distribution of positions inside a sphere, representing parasitic charged defects around the target boron-vacancy defect.To this end, we sample the elements of the position vector from standard normal distribution and multiply with d = r 3 √ x random distance, where x is sampled from a uniform random normalised distribution and r is the simulation sphere radius.Then we define an underlying lattice of bulk hBN and project the positions to the nearest lattice sites.We optimise two parameters, the charge density and the ODMR contrast, to fit to the experimental spectrum.The optimisation is done by the least squares method using parallelized brute force method on a two dimensional parameter grid refined in three consecutive optimisation cycles.

C. hBN crystals
The three hBN crystals used in this work were synthesised through the metal flux growth method described in Ref. [57], and irradiated at the Ohio State University Research Reactor, which produces a thermal neutron flux of 10 12 cm −2 s −1 .The crystal S1 is isotopically purified with 11 B (99.4%) and was irradiated with a neutron dose of 2.6 × 10 17 cm −2 .The two other crystals, S2 and S3, are purified with 10 B (99.2%) and irradiated with a dose of 2.6 × 10 16 cm −2 and 2.6 × 10 17 cm −2 , respectively.Neutron irradiation creates V − B centres (i) through damages induced by neutron scattering across the crystal and (ii) via neutron absorption leading to nuclear transmutation doping.The latter process strongly depends on the isotopic content of the hBN crystal since the neutron capture cross-section of 10 B is orders of magnitude larger than that of 11 B. As a result, neutron irradiation creates more efficiently V − B centres in hBN crystals isotopically enriched with 10 B.

D. Experimental details
The optical properties of hBN crystals are studied with a confocal microscope operating at ambient conditions.A laser excitation at 532 nm is focused onto the sample with a high numerical aperture microscope objective (NA=0.95).The PL signal is collected by the same objective and directed to a silicon avalanche photodiode operating in the single-photon counting regime.ODMR spectra are recorded by monitoring the PL signal while sweeping the frequency of a microwave field applied through a copper microwire deposited on the hBN crystal surface.

DECLARATIONS
• Data Availability The authors declare that the main data supporting the findings of this study are available within the paper and its Supplementary files.Part of source data is provided with this paper.The data that support the findings of this study are available from the corresponding author upon reasonable request.
• Code Availability The codes that were used in this study are available upon request to the corresponding author.Cao, Andrew Kauffman, and Kevin Herminghuysen for the irradiation services provided.We acknowledge the high-performance computational resources provided by KIF Ü (Governmental Agency for IT Development) institute of Hungary.
• Funding Open access funding provided by ELKH Wigner Research Centre for Physics.
• Author contribution PU carried out the DFT calculations and analysed the results under the supervision of AG.TCP, AD, BG, GC and VJ performed the experiments.JL and JHE provided the neutron-irradiated hBN crystals.PU, VJ and AG wrote the paper.All authors contributed to the discussion and commented on the manuscript.AG conceived and led the entire scientific project.
• Competing Interests The authors declare that there are no competing interests.

FIG. 1 .
FIG. 1. Geometry and electronic structure of the V − B centre in hBN. a Visualisation of the V − B centre and its ground state electronic spin density (yellow) in the mono-layer flake model used in the calculations.Green, blue, and white balls represent boron, nitrogen and hydrogen atoms, respectively.b Kohn-Sham level structure of the defect calculated with PBE0 functional.The orbitals are labelled according to the irreducible representations in D 3h point group symmetry.The subspace of occupied orbitals are indicated in the figure.As these lie close to each other, their labels are ordered in ascending energy from bottom to top.The spin density of the triplet ground state originates from the unoccupied doubledegenerate e orbital.

FIG. 2 .
FIG. 2. Response to an external strain component εxx. a Axial ZFS parameter D as a function of the external strain.bOrthorhombic ZFS parameter E as a function of the external strain.The cases of fixed and relaxed atomic positions are compared (blue and red points), as explained in the main text.c Local strain response for the externally applied strain.Blue dots show the one-to-one correspondence when the atomic positions are not allowed to relax.Red and green points show the enhanced local εxx strain and the activated local εyy strain around the defect, respectively, when its close vicinity is allowed to relax under the strain applied to the crystal.

FIG. 3 .
FIG. 3. Response to an external electric field component Ey. a Orthorhombic ZFS parameter E as a function of the external electric field.The cases of fixed and relaxed atomic positions are compared (blue and red points).b Estimated local strain from the piezo effect induced by the external electric field.

FIG. 4 .
FIG. 4. Simulated ZFS of the V − B centre in random perturbation environments.a ZFS distribution as a function of normal strain magnitude.b ZFS distribution as a function of electric field magnitude.
9 8 u 7 e o Z G 9 J e e S Z a z e g U n 1 5 B S h M b p I k o T x B W p 5 k 6 n m p n x f 7 m 3 d e e 6 m 5 X 9 H c z r 4 B

FIG. 5 .
FIG. 5. PL properties and ODMR spectra of hBN crystals with different concentration of V − B centres. a PL signal as a function of the optical excitation power for the three hBN crystals (S1, S2 and S3).The dashed lines are linear fit.b ODMR spectra recorded for S1, S2 and S3 crystals from top to bottom.The solid lines are the results of the fitting procedure based on the microscopic charge model described in the main text, from which we extract the charge density ρc.c Charge density ρc as a function of the relative density of V − B centres in the hBN crystal.The dashed line is a linear fit.

•
Acknowledgement This work was supported by the National Excellence Program for the project of Quantum-coherent materials (NKFIH Grant No. KKP129866), the Ministry of Culture and Innovation and the National Research, Development and Innovation Office within the Quantum Information National Laboratory of Hungary (Grant No. 2022-2.1.1-NL-2022-00004),the French Agence Nationale de la Recherche under the program ESR/EquipEx+ (grant number ANR-21-ESRE-0025), the Institute for Quantum Technologies in Occitanie through the project BONIQs and Qfoil.Support for hBN crystal growth is provided by the Office of Naval Research, awards numbers N00014-22-1-2582 and N00014-20-1-2474.The neutron irradiation was supported by the U.S. Department of Energy, Office of Nuclear Energy under DOE Idaho Operations Office Contract DE-AC07-051D14517 as part of a Nuclear Science User Facilities experiment.We acknowledge the support of The Ohio State University Nuclear Reactor Laboratory and the assistance of Susan M. White, Lei Raymond