Radiative properties of quantum emitters in boron nitride from excited state calculations and Bayesian analysis

Point defects in hexagonal boron nitride (hBN) have attracted growing attention as bright single-photon emitters. However, understanding of their atomic structure and radiative properties remains incomplete. Here we study the excited states and radiative lifetimes of over 20 native defects and carbon or oxygen impurities in hBN using ab initio density functional theory and GW plus Bethe-Salpeter equation calculations, generating a large data set of their emission energy, polarization and lifetime. We find a wide variability across quantum emitters, with exciton energies ranging from 0.3 to 4 eV and radiative lifetimes from ns to ms for different defect structures. Through a Bayesian statistical analysis, we identify various high-likelihood charge-neutral defect emitters, among which the native VNNB defect is predicted to possess emission energy and radiative lifetime in agreement with experiments. Our work advances the microscopic understanding of hBN single-photon emitters and introduces a computational framework to characterize and identify quantum emitters in 2D materials.


INTRODUCTION
Two-dimensional (2D) materials are rapidly growing as a platform for photon based quantum information technologies 1,2 . The discovery of single-photon emitters (SPEs) at point defects in hBN 3 has spurred an intense search for optically active defects in 2D crystals 4 . Compared to defects in bulk crystals, such as diamond and silicon carbide [5][6][7][8] , defects embedded in 2D materials promise to be more easily addressed and controlled. In the case of hBN, defect emitters exhibit a range of desirable properties, including high emission rate, room temperature stability, strong zero-phonon line (ZPL) and easy integration with other optical components using 2D hBN crystals 1,[9][10][11][12] .
A pressing challenge for defect emitters in hBN is identifying their atomic structure. Various possible structures have been proposed on the basis of density functional theory (DFT) calculations and their comparison with experiments 3,[13][14][15][16][17][18][19] . However, while DFT can provide valuable insight into the formation energy, symmetry and electronic structure of SPE defects, it cannot address key aspects of point-defect SPEs such as their excited states and radiative processes responsible for light emission. To complicate the matter further, similar to other 2D materials and their defects [20][21][22] , optical transitions at defects in hBN are dominated by excitonic effects 23 , which require specialized first-principles calculations beyond the scope of DFT. Quantum chemistry approaches and density matrix renormalization group have also been used to investigate specific defect structures 24,25 , but wide comparisons among different structures are still missing.
The Bethe-Salpeter equation (BSE) can accurately predict optical properties and excitons in materials 26 . It also enables, due to recent advances, precise calculations of radiative lifetimes in 2D and bulk crystals [27][28][29][30] . The radiative lifetime plays an important role in the study of SPEs as it determines the shortest decay time constant in the second-order photon correlation function 31,32 and can also be measured directly from fluorescence intensity decay 3,10 . Applying the BSE and related methods to defect emitters in hBN would enable direct comparisons between theory and experiment of the emission energy and radiative lifetime, providing valuable information to identify defect SPEs. Yet, firstprinciples BSE calculations are computationally costly to carry out on defect structures, and the radiative lifetime calculations are only a recent development.
In this work, we employ the BSE approach to compute from first principles the optical properties, transition dipoles, excitons and radiative lifetimes of atomic defects in hBN. We examine a large pool of candidate SPE structures, spanning native defects and carbon or oxygen impurities, to correlate their atomic structures with their photophysics. We find that different quantum emitters exhibit radiative lifetimes spanning six orders of magnitude and emission energies from infrared to ultraviolet. Bayesian statistical analysis is employed to correlate our results with experiments and identify high-likelihood charge-neutral SPEs in hBN, among which we find the V N N B defect to have the highest likelihood. In-depth calculations on the V N N B defect highlight the strong dependence of its radiative properties on small perturbations to its atomic structure. The dependence of the defect radiative properties on dielectric screening is analyzed by comparing monolayer and bulk hBN results. Our systematic investigation addresses key challenges for characterizing the excited states and radiative properties of defect emitters in 2D materials. calculations 33 , followed by BSE calculations to obtain the exciton energies and wave functions and from them the optical absorption, transition dipoles and radiative lifetimes (see the "Methods" section). In the following, we denote the defects in hBN as X N Y B if neighboring N and B atoms are replaced by species X and Y, respectively, where X and Y can be a vacancy (denoted as V) or another element 3,13 . We focus on emitters in the interior of the 2D crystal 34,35 and do not consider defects that would likely appear at the sample edges or corners 36,37 .
The electronic energies obtained using DFT, while in general not representative of electronic or optical transitions, can be used for estimating qualitative trends. Figure 1 shows the lowest (HOMO-LUMO) spin-conserving transition energy of the candidate defects, obtained from DFT, together with the emission polarization inferred from structural symmetry. The defect structures considered here exhibit three different types of local symmetries, D 3h , C 2v , and C s . In the high symmetry D 3h structure, adopted by N B , V N , B N , C B , C N , O N , emitted light cannot be linearly polarized. Conversely, linearly polarized emitted light, as observed experimentally in hBN SPEs 11 , is possible in the C 2v and C s symmetries. In the C 2v configuration, which is the most common among the defects investigated here, the threefold rotational symmetry is broken but all the atoms remain in-plane, preserving the mirror symmetry with respect to the crystal plane. The C s symmetry found in the V N N B , V N C B , V N O B and O B defects is instead associated with an out-of-plane distortion that breaks the mirror symmetry about the plane. The DFT transition energies for the 22 defect structures range from 0 to 3.5 eV. In contrast, the ZPL of the measured SPEs are in the 1.6−2.2 eV energy range 11 , as shown by the shaded region in Fig. 1. While candidate structures with D 3h symmetry can be ruled out, C 2v and C s structures with exceedingly small or large DFT transition energies also appear unlikely on the basis of the DFT results.
Starting from the DFT ground state, for selected defects we compute the excited state properties with the GW-BSE method, obtaining the quasiparticle energies in the one-shot G 0 W 0 approximation and the exciton energies and wave functions with the BSE, which captures electron-hole interaction and excitonic effects. We combine the solutions of the BSE with an approach we recently developed 29,30 to compute the radiative lifetime of an exciton state from Fermi's golden rule. Generalizing our previous formula for isolated (0D) emitters 29 to include anisotropic dielectric screening in hBN, the radiative decay rate γ S of an exciton state S is: where ϵ xy (k xy ) and ϵ z are the in-plane and out-of-plane dielectric function of hBN, respectively, k xy is the in-plane photon wavevector, E S is the exciton energy and p S,xy and p S,z are the corresponding components of the exciton transition dipole. The radiative lifetime τ S below is defined as the inverse of γ S and accounts only for the radiative decay channel. We use a constant in-plane dielectric function for bulk hBN (ϵ xy = 5) at optical wavelengths, and for monolayer hBN we take into account the dependence on wavevector k as ϵ xy (k) ≈ 1 + 2πα 2D |k| 38 , where α 2D is a constant equal to 0.4 nm 39 . In this approach, which is appropriate for 2D materials, the in-plane dielectric function of monolayer hBN reduces to a value of 1 at optical photon wavelengths. Figure 2 shows the computed radiative lifetimes and lowest bright exciton energies for nine selected defects, including V N , B N , We find that the exciton energy can differ significantly−by as much as 1 eV− from the corresponding DFT transition energy, which fails to account for screening and electron-hole interaction effects. This result is a testament to the importance of excited state calculations for investigating SPEs. Surprisingly, we find that the computed  Writing the exciton transition dipole in Eq. (1) as jto highlight its physical meaning of an atomic-scale dipole, and setting r S j j equal the B-N bond length, we obtain an upper bound estimate of the exciton transition dipole (of 6.9 Debye) because the excitons considered here are tightly bound to the defect. The resulting radiative lifetime as a function of exciton energy, as shown by the blue line in Fig. 2, gives an approximate lower bound to the calculated radiative lifetimes. Our calculated exciton energies and radiative lifetimes clearly follow the trend set by this lower bound. The physical insight of this analysis is that due to incomplete overlap of the electron and hole wavefunctions, the exciton transition dipole for most defects is significantly smaller than the bond length, leading to longer radiative lifetimes than this theoretical bound. We note that the experimental emission energies and lifetimes for defect emitters in hBN, indicated by the gray area in Fig. 2, fall almost entirely below this lower bound, suggesting that nonradiative decay processes need to play an important role if atomic point defects like the ones discussed here are responsible for SPE.
Bayesian analysis of the likelihood of candidate defect structures As the candidate defect structures have properties distributed across a wide range, it is challenging to pinpoint the correct defect structure from our results. As such, we pursue a quantitative analysis of the relative likelihood of the various structures using Bayesian inference, a statistical approach for dealing with uncertainties and combining information from different categories, in which the likelihood of a hypothesis is updated as more evidence or information becomes available 40 . The workflow of our analysis is shown in Fig. 3. We first generate candidate structures by systematically enumerating all possible structures where n adjacent atomic sites are replaced by a vacancy, antisite or impurity point defect. We then obtain the list of defects studied in our work by considering structures with n = 1 or 2 replacements, as defects spanning multiple sites are typically less likely to occur. We use carbon and oxygen as the impurity elements, the most likely impurities in hBN according to the current literature. This list of candidates includes structures similar to known solid-state quantum emitters such as the NV and SiV centers in diamond, and includes most of the hBN defect emitters proposed in the literature so far [13][14][15][16][17][18][19] as well as additional structures that were not previously considered.
A prior likelihood p(h) is then assigned to each candidate structure h as an initial guess of how likely it is that the defect structure can appear in the hBN crystal. The prior likelihood can be tailored to describe various experimental scenarios. For example, if carbon impurities are added in the experiment, as shown in recent work 41 , one could include the presence of a carbon impurity in the prior likelihood, therefore favoring the posterior likelihood of carboncontaining emitters. Here, without referring to specific experimental conditions, we choose the prior likelihood p(h) on the basis of the structural complexity of the defect, using the single exponential form pðhÞ / S h A ÀC h , where C h measures the complexity of the structure, defined as the number of vacancies and antisite defects plus twice the number of impurities, A is a constant, and S h is a degeneracy factor due to symmetry-related configurations.
We subsequently combine our computational results with experimental evidence to obtain a posterior likelihood p(h|E) for each structure 40 : where the likelihood function p(E|h) is the likelihood that the structure h is compatible with a specific piece of experimental evidence E on the basis of our calculations. The resulting posterior likelihood, which is updated in multiple steps using several computed properties, quantifies which structures among those considered are more likely to be defect emitters in hBN. The properties of the candidate SPE defects and the results of our Bayesian inference analysis are summarized in Table 1. Moving in the table from left to right, the prior likelihood p(h) is updated to posterior likelihood p(h|E) using calculations presented in subsequent columns. The first computed posterior likelihood, p 1 (h|E), uses the DFT results to infer whether the defect symmetry is compatible with linearly polarized emission and whether the DFT transition energy, within a standard deviation of~0.5 eV, falls in the 1.6−2.2 eV emission range found in experiments 10,11,34 . The final posterior likelihood p 2 (h|E) in the rightmost column of Table 1 gives the likelihood of the most likely defect structures when all factors have been taken into account, including comparison between the lowest bright exciton energy and the experimental  Fig. 3 Bayesian inference workflow. Candidate structures are initially generated through combinatorial enumeration of point defects in hBN and attributed a prior likelihood based on structural complexity. The likelihood of each defect is then updated using results from our DFT calculations. For the most promising defect structures, the likelihood is refined using GW-BSE calculations.

Posterior Likelihood p(h) p1(h|E) p2(h|E)
emission energy, and comparison with experiment of the radiative lifetime from our GW-BSE calculations. We emphasize that the absolute value of the likelihood is arbitrary and a value of 1 does not imply certainty. The likelihood shown in Table 1 is normalized so that the maximum value on each column is 1. The detailed rationale of our Bayesian analysis is described in the Supplementary Note 3. The main finding from the data in Table 1 is that the V N N B defect, which was originally proposed as a SPE in hBN on the basis of DFT calculations 3,13,18 , possesses optical and radiative properties that most accurately match the experimental results, even after taking into account excitonic effects and radiative lifetimes. The next most likely structure is the oxygen impurity defect, O B , with an emission energy of 1.52 eV, just below the experimental range. On the other hand, we find that the V N C B defect, which is also considered a likely candidate in the literature based on DFT calculations 13,18 , has an emission energy of 2.6 eV that lies too far above the experimental energy range when excitonic effects are included, making it a less likely candidate. Although our analysis excludes V N C B and other defects with higher emission energies as candidates for the 1.6−2.2 eV SPEs, these structures may still be good candidates for SPEs in the ultraviolet range 42,43 . As the most likely defect emitter in hBN is the V N N B structure according to our analysis, we investigate it in more detail to better understand its radiative properties. Figure 4(a) shows that the relaxed atomic structure of the V N N B defect in monolayer hBN exhibits an out-of-plane displacement of the central nitrogen (N) atom by z = 0.66 Å. If the structure is made planar by constraining the N atom to the hBN plane, the total energy of the system is only moderately higher (by 0.125 eV) than the equilibrium structure. Therefore the total energy forms a shallow double-well potential as a function of the out-of-plane N atom displacement. Consistent with recent findings 44 , this small energy difference suggests that the structure is soft in the out-ofplane direction and could fluctuate due to external perturbations.
To take this result into account, we compute the optical absorption spectrum of the V N N B defect for out-of-plane displacements z of the N atom in the 0−0.83 Å range, interpolating between the in-plane and out-of-plane equilibrium structures to obtain the intermediate metastable structures. We find that the absorption spectrum changes drastically even for small changes of z, as shown in Fig. 4(c). The spectrum is dominated by two exciton absorption peaks at low energy. At the equilibrium position (z = 0.66 Å), the first peak at 1.92 eV is associated with a transition from the top valence to the bottom conduction spin-majority electronic states, whose wave functions are shown in Fig. 4(b). As the displacement z of the N atom is decreased from the equilibrium value to zero, which corresponds to a planar structure, the energy of the first exciton peak decreases at first but then increases again at small z values while gaining oscillator strength monotonically. The second peak at 2.7 eV is also due to a transition from the top valence to the bottom conduction spin-minority bands with relatively large oscillator strength at the out-of-plane equilibrium N atom position. Both the energy and the oscillator strength of this peak decrease monotonically as z decreases [see Fig. 4(c)]. For the planar structure (z=0), the second peak becomes the lowest-energy transition but completely loses its oscillator strength and becomes a dark state. These significant changes in the two lowest-energy exciton peaks are accompanied by substantial changes in the exciton radiative lifetimes. Although the lifetime of the lowest exciton is 334 ns for the DFT equilibrium out-of-plane distortion, it becomes only 19 ns if the V N N B structure is kept planar (z = 0 case), as shown in Fig. 4(d). This value of the radiative lifetime for the planar structure is within a factor of 2−3 of experimental results 11 . Because many experiments are carried out on thicker hBN samples, we also calculated the radiative lifetime of the V N N B defect embedded in bulk hBN (instead of monolayer), a setting resembling more closely SPE experimental measurements in thicker hBN samples (see the "Methods" section). Due to interlayer interactions, the out-of-plane displacement of the N atom is reduced to 0.45 Å in bulk hBN. Going from monolayer to bulk reduces the radiative lifetime of the V N N B defect to 147 ns due to changes in the dielectric environment and transition dipole, both of which affect the radiative emission rate in Eq. (1). This result shows that different dielectric environments, for example due to the use of different substrates, may have a strong impact on the radiative properties of the defect. Therefore we conclude that the radiative lifetime of the V N N B defect is highly sensitive to both outof-plane structural distortions and the dielectric environment.

DISCUSSION
While our computed lifetimes are in the 20−350 ns range for different V N N B structures, these values pertain to an ideal defect with 100% quantum yield, whereas defect emission measured so far in hBN is likely associated with a lower quantum yield. To explain the discrepancy between the experimental and calculated lifetimes, note that nonradiative processes, which are not included in our calculation of intrinsic radiative lifetimes, could significantly reduce the apparent lifetime measured in experiments. For example, a defect with a measured 10 ns lifetime but only a 10% quantum yield would correspond to a 100 ns intrinsic radiative lifetime in the same range as our computed values. Future SPE measurements in hBN taking into account the quantum yield may provide a better estimate of the intrinsic SPE radiative lifetime and enable a more accurate comparison with our calculations. In our study, we did not take into account the Frank-Condon shift or exciton-phonon coupling effects such as zero-point renormalization, which are beyond the scope of this work. Although these effects may cause small shifts of the exciton energy and lead to broadening or satellite peaks 45,46 , the overall transition dipole remains the same within the Frank-Condon approximation, and thus we expect minimal changes in the radiative lifetimes due to these effects. Note also the possible presence of charged defects in experiments 4 , which are not considered in this work.
In summary, we have investigated the excited state and radiative properties of many candidate defect SPEs in hBN. Our calculations address the photophysics of these defect structures, including their optical transitions and radiative lifetimes, using calculations that accurately account for excitonic effects and anisotropic dielectric screening. Our Bayesian statistical analysis allows us to identify the native V N N B defect as the most likely candidate within a wide pool of over 20 charge-neutral native defects and carbon or oxygen impurities. This work explores the fertile ground at the intersection of excited-state first-principles calculations and Bayesian learning methods, and the application of this framework to quantum technologies. The search for solidstate quantum emitters will greatly benefit from such precise firstprinciples calculations of excited states in materials.

Details of the first-principles calculations
We carry out DFT calculations in the generalized gradient approximation using the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional 47 . Our spin-polarized DFT calculations employ the plane-wave pseudopotential method implemented in QUANTUM ESPRESSO 48 . The defects are placed in a 5 × 5 supercell of hBN with the lattice constant kept fixed at 5 × 2.504 Å while the atoms are fully relaxed without symmetry constraints. For calculations on monolayer hBN, a vacuum of 15 Å is used to avoid spurious inter-layer interactions with the periodic replicas. We use ONCV pseudopotentials 49,50 for all atoms, a plane-wave kinetic energy cutoff of Fig. 4 Properties of the V N N B defect. a Atomic structure in top and side views, with N atoms shown as blue and B atoms as green spheres. Also shown in the side view is the N atom out-of-plane displacement z. b GW electronic band structure computed using a 5 × 5 supercell. The wave function of the two majority-spin bands that mainly contribute to the lowest-energy exciton are plotted on the right. c Absorption spectrum and (d) radiative lifetime as a function of the out-of-plane displacement z of the central N atom. 80 Ry and a 3 × 3 k-point Brillouin zone grid. The GW-BSE calculations are carried out with the Yambo code 51,52 using a 2D slab cutoff of the Coulomb interaction. For calculations of defects in monolayer hBN, a 3 × 3 k-point grid is employed together with an energy cutoff of 10 Ry for the dielectric matrix. The number of empty bands included in the GW calculation (for the polarizability and self-energy summations) is 7 times the number of occupied bands. At least 6 valence and 6 conduction bands are included when diagonalizing the BSE matrix. These parameters converge the exciton energy to within 0.1 eV in our test calculations, which are shown in the Supplementary Note 1. For calculations on defects in bulk hBN, we use a 5 × 5 × 1 supercell containing two AA'-stacked layers. The bulk case uses a 2 × 2 × 3 k-point grid; all other parameters are the same as in the monolayer case.

DATA AVAILABILITY
All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Information. Additional data related to this paper may be requested from the authors.