Vibrational signatures for the identification of single-photon emitters in hexagonal boron nitride

Color centers in h-BN are among the brightest emission centers known yet the origins of these emission centers are not well understood. Here, using first-principles calculations in combination with the generating function method, we systematically elucidate the coupling of specific defects to the vibrational degrees of freedom. We show that the lineshape of many defects exhibits strong coupling to high frequency phonon modes and that C$_{\text{N}}$, C$_{\text{B}}$, C$_{\text{B}}$-C$_{\text{N}}$ dimer and V$_{\text{B}}$ can be associated with experimental lineshapes. Our detailed theoretical study serves as a guide to identify optically active defects in h-BN that can suit specific applications in photonic-based quantum technologies, such as single photon emitters, hybrid spin-photon interfaces, or spin-mechanics interfaces.


I. INTRODUCTION
Hexagonal boron nitride (h-BN) is a wide band gap van-der-Waals solid. In its exfoliated form, h-BN is a stable two-dimensional material that retains its wide band gap of about 6 eV. This large band gap supports a diversity of optically active defect centers, exhibiting a wide range of emission energies between 1.2 and 5.3 eV [1,2]. Recently, it was shown that mono and few-layer h-BN can host room temperature single-photon emitters (SPEs) [2], which sparked enormous interest in the field of photonic-based quantum technologies [3,4]. Follow-up experiments studied SPEs in h-BN in more detail [2,[5][6][7][8][9][10] and, remarkably, such emitters were shown to exhibit Fourier transform-limited emission up to room temperature [11]. An understanding of the optical properties of SPEs in h-BN would enable selecting specific defect centers to serve as single photon sources [4], as hybrid spin-photon [12] or spin-mechanics interfaces [13,14].
In order to take full advantage of these SPEs, it is paramount to identify the responsible defect structures and corresponding electronic transitions. However, so far no consensus has been reached concerning the origin of single-photon emission. While the spatial localization of the single-photon emission strongly suggests point defects to be the culprit, the assignment to a specific defect in h-BN is hampered by the range of zero-phonon lines (ZPLs) that are observed and the large number of potential defect candidates.
For ZPLs between 1.6 and 2.3 eV, in most cases pronounced phonon sidebands (PSBs) at about 165 meV below the ZPL are observed [2,[5][6][7][8][9]. In some experiments, PSBs for some SPEs exhibit, however, a double-peak structure at around 160 meV and 190-200 meV [9,10]. The similarity of the PSBs across measurements suggests the emission to be due to either multiple different defects with similar geometry or a single defect with variable excitation energy. In a recent experiment, the variable excitation energy was attributed to strain effects [10], but Stark shifts have been suggested as well [15]. It has also been proposed that there are two families of emitters around 2 eV with different electronic structure and ZPLs of 1.88 and 2.14 eV, respectively, that can be distinguished not only based on their ZPL but also their quantum efficiency [16].
h-BN also exhibits luminescence in the ultraviolet [17], with a recent experiment [18] demonstrating single photon emission at 4.1 eV. The structure of this defect remains unknown but has been proposed to originate from carbon defects, although the emission intensity does not exhibit correlation with C impurity concentration [19]. The highest intensity peak of the PSB of the 4.1 eV emitter has been suggested to originate from coupling to a 187 meV local phonon mode [17] and coupling to the zone center longitudinal optical (LO) phonon mode at 200 meV [20].
Recent theoretical studies focused on calculating the electronic structure of point defects in monolayer and bulk h-BN based on first-principle methods [21,22]. The ZPLs for rather many defects have been calculated within the accuracy permitted by current density functional theory (DFT) methods based on hybrid functionals. Furthermore, the PSBs have been analyzed using phenomenological models [9,23] using a few selected phonon modes. To the best of our knowledge, a combined defect and PSBs study has, however, only been performed on the defects N B -V N and C B -V N [21,24]. In the latter studies it was concluded that the calculated emission spectra of N B -V N do not agree with the measured lineshapes but that C B -V N might be a SPE. Besides these two defects, information about how specific defects couple to the vibrational degrees of freedom is scarce.
In the present work, we contribute to closing this knowledge gap concerning point defect-related emissions in h-BN by considering both charged and charge neutral transitions and the resulting emission lineshape for a set of the most common intrinsic and extrinsic point defects. Importantly, we calculate the combined defect and PSB emission spectrum and, thus, can assess the vibrational fine structure. The resulting emission spectrum can be used as an experimentally accessible fingerprint to identify defect-related SPEs. This is possible since the vibrational fine structure of the emission spectrum due to an electronic transition on a point defect is highly sensitive to changes in the local distortion between initial and final state [25]. In the following, we consider vacancies (V N , V B ) and antisites (N B , B N ) as the most important intrinsic defects as well as carbon impurities (C B , C N ), vacancy-impurity complexes (C B -V N , C N -V B ) and one antisite-vacancy complex (N B -V N ). This selection covers most of the defects that have been proposed as SPEs in h-BN.

A. Defect formation energy
The formation energy of a defect in charge state q is given by where E defect and E ideal are the total energies of the defective and ideal systems, respectively. ∆N i denotes the change in the number of atoms of type i between defective and ideal system, while ε VBM and ∆µ e are the valence band maximum (VBM) position and the (relative) electron chemical potential, respectively. The chemical potentials µ i of B and N are coupled to each other via where E BN is the cohesive energy of h-BN. [26] Below, we consider the nitrogen-rich limit, where µ N = − 1 2 E N2 and the boron-rich limit, where µ B is taken as the negative cohesive energy of elemental α-B (spacegroup R3m). The chemical potential of carbon is set to the one of graphene throughout. The charge transition level (CTL) between charge states q and q is the value of the electron chemical potential for which the formation energies of the defect in charge state q and q are equal. Throughout this study, CTLs are reported with respect to the VBM.

B. Lineshape of emission spectrum
PSBs arise due to emission from the electronic excited state to the vibrationally excited electronic ground state. The structure of the emission spectrum can be computed from a knowledge of the phonon spectrum and the difference in the ionic configurations associated with excited and ground states ∆R. The extent of the lattice distortion can be measured by the magnitude of the mass weighted difference in ionic displacements where the sum runs over all atoms in the defect cell. In this work, the lineshape is computed using the generating function approach [25,27,28]. The central quantity in the generating function approach is the electronphonon spectral function, which depends on the coupling between lattice displacement and vibrational degrees of freedom. The latter information is encoded in the so called partial Huang-Rhys (HR) factor This expression is obtained in the low temperature and parallel mode approximation, i.e., the frequencies and eigenvectors of both ground and excited electronic states are related by a simple translation. Q k is the projection of the lattice displacement on the normalized collective displacement described by phonon mode k and given by [29] Q k = a √ m a n a k ∆R a , where a runs over all the atoms in the computational cell and n k is the normalized ionic displacement vector corresponding to phonon mode k.
The electron-phonon spectral function is then obtained by summation over all modes It has dimensions of inverse energy in the same units as ω. The integral over the electron-phonon spectral function gives the (total) HR factor of the transition. The electron-phonon spectral function is then transformed to the time domain to obtain S(t) = dωS(ω) exp(−iωt) and the generating function G(t) is obtained by exponentiation of S(t) Fourier transformation of the generating function yields the lineshape function Here, κ is a broadening parameter that governs the width of the ZPL. It is necessary for numerical reasons and does not represent thermal broadening. It should be chosen as small as possible to minimize its effect on the spectrum and as large as necessary top achieve a smooth representation of the spectral functions. Its role is thus akin to the smearing width that is adopted when computing electronic and phonon densities of states. Similarly, there is a correlation between the convergence with respect to smearing and the number of modes (i.e. the size of the supercell in the present case or the density of the Brillouin zone mesh when computing a density of states). Here, we use a value of 6 meV, which is chosen to achieve the balance described above (see Fig. S11 in the Supplementary Information (SI)). Finally, the lineshape function is related to the luminescence intensity via where C is a normalization constant chosen such that dω L(ω) = 1. The localization of a phonon mode can be measured using the inverse participation ratio (IPR) [29] which can assume values between 1 and N , where N is the number of atoms in the computational cell for which the phonon spectrum has been calculated. Smaller and larger values indicate more and less localized character.

C. Computational details
All DFT calculations were carried out using planewave basis sets [30] and the projector augmented wave method [31,32] as implemented in the Vienna ab initio simulation package [33] (VASP, version 5.4.4). Exchangecorrelation contributions were obtained using the semilocal PBE functional [34] and the hybrid HSE06 functional [35] using both the standard value for the mixing parameter α = 0.25 as well as a value of α = 0.60 tuned to reproduce the band gap as detailed below (Sect. III A). A plane wave basis set with a cutoff energy of 550 eV was employed to represent the electronic wave functions. Geometry optimization was performed for all systems, during which the atomic positions were allowed to relax until all forces were less than 20 meV/Å. Brillouin zone sampling was performed using a grid of 21 × 21 × 1 for the primitive hexagonal (2-atom) unit cell. Defect calculations were carried out using a 8×8×1 supercell with 19.88Å vacuum while the Brillouin zone was sampled using the Γ-point only.
The ZPLs arising from transitions between defectinduced levels and band states were computed using the HSE06 hybrid functional using a mixing parameter of α = 0.6 to correct for the band gap error of semi-local functionals (see Sect. III A for the motivation of this choice of mixing parameter). For charged defect cells, a correction of q 2 0.5 eV was added to the total energy to account for periodic image charge effects and potential offsets [36,37]. This approach yields CTLs that agree well (within 30 meV) with the values presented in the erratum of Ref. 36.
Charge neutral excitations were modelled using the ∆SCF method, in which the electronic occupations are constrained. The description of electronic states that are localized at the level of semi-local exchange-correlation functionals such as PBE is commonly only weakly affected by the addition of exact exchange. (For illustration, PBE and standard HSE06 (α = 0.25) yield ZPLs for V −1 B of 1.71 eV and 1.84 eV, respectively.) Since in the present case the convergence of excited state calculations with hybrid functionals is both very tedious and computationally demanding, all charge neutral defectdefect transitions were treated at the semi-local DFT level (PBE).
Vibrational spectra were obtained at the semi-local DFT level (PBE) using the phonopy package [38]. In the computation of the electron phonon spectral function S(ω), the Kronecker δ in Eq. (6) was approximated using normalized Gaussians with a smearing of 6 meV. The integration over the time domain in the Fourier transform to obtain A(ω) was performed over 2 ps with a time spacing of 1 fs. Convergence of the spectral distribution function is demonstrated in the SI.

A. Pristine h-BN
Without taking into account zero-point effects, the PBE functional yields a lattice parameter of 2.51Å for h-BN while one obtains 2.49Å with HSE06 using α = 0.25. These values compare well with the experimental lattice parameter for bulk BN of 2.51Å at 10 K [39]. For consistency, all calculations below, including those based on the standard HSE06 functional (mixing parameter α = 0.25), were carried out using a lattice constant of 2.51Å.
The calculated electronic (or single-particle) band gap, on the other hand, measures 4.67 and 5.71 eV with PBE and standard HSE06 (α = 0.25), respectively, in line with other theoretical studies [40]. The experimentally measured optical band gap is 6.1 eV and 6.0 eV for monolayer h-BN [41] and bulk BN [42], respectively. While these numbers appear rather close to the value obtained using standard HSE (α = 0.25), a sound comparison requires accounting for the exciton binding energy, which is very large both in bulk [43] and monolayer h-BN [40], as well as the effect of zero-point motion via electron-phonon coupling [44]. It is therefore more instructive to compare the DFT values with the results from G 0 W 0 calculations [40], which provide a value of 7.1 eV, corresponding to the single-particle band gap in the absence of both exciton formation and zero-point motion. To reach this value using the HSE06 functional one requires a mixing parameter of α = 0.60, yielding a band gap of 7.2 eV. Below, we therefore use HSE06 with a mixing parameter of α = 0.60 to compute defect formation energies and ZPLs derived thereof.
Before considering the vibrational spectra of defects, a closer inspection of the vibrational spectrum of pristine h-BN is instructive. The vibrational spectra from PBE and HSE06 are very similar, with the latter yielding a slightly stiffer response overall (Fig. 1). [45] Given the 2D character of h-BN, the phonon density of states (DOS) can be decomposed into an in-plane and out-ofplane part. In part due to the quadratic dispersion inherent to 2D materials [46], the lower frequency part of the spectrum is dominated by out-of-plane modes. The out-of-plane partial DOS features two pronounced bands ranging from 0 to about 40 meV and from approximately 65 to 100 meV (values from PBE), respectively, with two pronounced peaks at 40 and 80 meV. The in-plane partial DOS covers the entire frequency range spanning up to 187 meV in good agreement with other first-principle studies [47,48]. The most notable feature is an asymmetric peak at 154 meV.

B. Defect energetics
From the outset we considered vacancies (V N , V B ) and antisites (N B , B N ) as possibly relevant intrinsic defects as well as carbon impurities (C B , C N , C B -C N ), vacancyimpurity complexes (C B -V N , C N -V B ) and one antisitevacancy complex (N B -V N ), see Fig. 2a for atomic structures. To determine the energetically most favorable defects and defect configurations along with pertinent charge states we computed the defect formation energies under both B and N-rich conditions (Fig. 2b). Where comparable the results are consistent with previous work on monolayer h-BN [36]. Under N-rich conditions N B (donor with q = +1, 0) and V B (acceptor with q = 0, −1) are the most favorable intrinsic defects, whereas under B-rich conditions B N (acceptor, q = 0, −1) and V N (ambipolar with q = +1, 0, −1) have the lowest formation energies. With regard to extrinsic defects, we find C B (q = +1, 0), C N (q = 0, −1), and C B -C N (q = 0, −1) to be lowest energy defects under both B and Nrich conditions. The defect complexes involving substitutional impurities and vacancies, specifically C B -V N and N B -V N , have high formation energies but are nonetheless included in the analysis below since they have been discussed as potential SPEs before [2,21,49].

C. Transition types
In the following, we consider two types of transitions: (i) localized-to-delocalized (LD) transitions involve a localized defect level and a (delocalized) band edge state; (ii) localized-to-localized (LL) transitions involve two defect levels both of which reside inside the band gap when they are occupied.
To illustrate the features and emphasize similarities and differences between these transition types, it is instructive to recall the relations between configuration coordinate diagram, defect formation energies, CTLs, and defect levels.
In the case of LD transitions, illustrated here by C N , the absorption (emission) energy corresponds to the 0/−1 CTL obtained when the atomic configuration is constrained to the equilibrium C 0 N (C − N ) geometry (Fig. 3a). The ZPL then corresponds to the CTL obtained for the equilibrium geometries in either charge state (Fig. 3b). If the electron chemical potential ∆µ e coincides with the VBM, emission occurs by capturing a hole from the VB edge, depleting the defect level associated with C N (Fig. 3c). If the electron chemical potential resides at the conduction band minimum (CBM) the initial and final states are inverted. Since a delocalized band state is involved in this transition type, a proper description of the position of the band edges is crucial.
To illustrate LL transitions, we consider here the C B -C N defect, which presents a particular simple two-level system. Emission occurs from the excited state C B -C 0, * N , in which the highest occupied defect state resides above the lowest unoccupied level, to the ground state C B -C 0 N (Fig. 3d). The difference in character of the highest occupied defect level in ground and excited states implies a difference in local potential that gives rise to a considerable lattice relaxation after emission (Fig. 3e) that underlies the Stokes shift.
The possible LD and LL transitions involving the defects considered here are compiled in Table I.

D. Transitions on intrinsic defects
Next we turn to a survey of the transition energies and subsequently an assessment of the lineshapes of intrinsic defects. The prominent high-frequency PSB of emitters in the 1.6 to 2.3 eV range and at 4.1 eV indicates an effective phonon frequency of ω eff ∼ 100 meV. The effective frequency is coupled to the HR factor S and the lattice distortion connecting the initial and final configurations ∆Q according to  As experimentally measured HR factors fall in the range of S = 1 to 2, the lattice distortion of potential emitters must therefore be relatively small. For example, a transition coupled to an effective mode with a frequency of 100 meV must have ∆Q <0.41 √ amuÅ in order to have an HR factor below two.
Most transitions on intrinsic defects exhibit, however, rather large lattice distortions with ∆Q well above 1 √ amuÅ (see ∆Q in Table I) that disqualify them as possible narrow band emitters. These defects often exhibit geometries in which one or several atoms are located outside of the BN plane, which can be related to the large structural differences between different electronic states. Transitions involving V N and N B are therefore not considered further.
B N can host a LL transition with a ZPL of 1.14 eV according to PBE with a relatively small ∆Q. As the transition is rather far from the spectral range of interest here (> 1.6 eV), transitions on B N are not considered further either. N B -V N can be ruled out as well based on the large value of ∆Q. The one remaining intrinsic defect V B is considered in detail in the following section.

Boron vacancy (VB)
The symmetry of V 0 B is C 2v , while the symmetry of V − B is D 3h . The 0/ − 1 CTL resides 4.12 eV above the VBM, which is much larger than an earlier theoretical value of 2.4 eV [50]. The latter though was obtained using a much smaller mixing parameter of α = 0.32, which as discussed in Sect. III A yields a misleading agreement with the experimentally measured optical band gap. The Here, we consider only the highest occupied molecular orbital (HOMO)/lowest unoccupied molecular orbital (LUMO) transition in both the D 3h and C 2v symmetries of the excited state with ZPLs of 1.72 eV and 1.63 eV, respectively. The difference of ∼90 meV between the considered ionic configurations of the excited state suggests that C 2v should be the equilibrium configuration. Below we, however, consider the lineshapes of both the configurations with D 3h and C 2v symmetry, due to the small energy difference and the possibility that semilocal DFT provides an incorrect description of the equilibrium ionic configuration of the excited state.
First, the case where the initial state symmetry is C 2v is considered with a total HR factor of 2.45. In this case there is a strong coupling of the electronic transition to a single mode at 26 meV with a partial HR factor of 0.66 (27% of the total HR factor) (the spectral function is shown in the SI, Fig. S9). The IPR of this mode is 92, suggesting that it is likely a delocalized mode (the maximal value of the IPR that can be reached in this supercell is 128). At energies around 40 meV, there are 3 modes with partial HR factors of 0.11 to 0.29 with IPRs of 40 to 50. In the high-frequency end, there is a coupling to a 162 meV mode with a partial HR factor of 0.05 (2% of the total HR factor). The resulting normalized emis-sion lineshape is shown in Fig. 8, where the spectrum has been broadened by tuning the damping parameter (κ) in Eq. 7 (see Fig. S11 for an illustration of how the spectral distribution function changes with increasing κ in the case of C N ).
Next, the emission lineshape of the transition from the D 3h symmetric initial state is considered (Fig. 4a), which has a much smaller total HR of 0.91. The spectral function (Fig. 4b) has two prominent peaks at 40 meV and 162 meV. The 40 meV peak results from the coupling to a single phonon mode and the coupling to other phonons with energy in the vicinity of 40 meV is very weak. The IPR of the 40 meV mode is 50 and the partial HR of this mode is 0.34 (37% of total HR). The phonon eigenvector for the 40 meV and the 162 meV mode is shown in Fig. 4b. The IPR of the 162 meV mode is 26 and the partial HR factor is 0.07. While this value is slightly larger than the corresponding value in the emission from the C 2v initial state, the relative contribution is much larger at 8% of the total HR factor. Therefore, the peak at 162 meV in the emission lineshape is much more pronounced in the case of emission from the D 3h state. The normalized emission lineshape of the emission from the D 3h is shown in Fig. 8.

E. Transitions on extrinsic defects
Compared with intrinsic defects among which only V B is a viable candidate for single-photon emission, extrinsic defects involving carbon feature a multitude of suitable electronic transitions (Table I). Transitions involving C B -V N (S) and C N -V B exhibit large lattice distortions and are therefore excluded from further analysis as argued above. In the following we therefore focus on C N , C B , C B -V N (T), as well as C B -C N dimers.  Fig. 3c). Arrows indicate the direction of a transition. Singlet and triplet states are marked by (S) and (T ), respectively. Localized-delocalized transitions can proceed in either direction, whereas localized-localized transitions always commence from the excited state, the latter being marked by asterisks. The "Modes" column identifies which phonon spectrum was used for the computation of the lineshape function. ZPL: zero-phonon line; ∆Q: magnitude of the lattice distortion between the two defect configurations involved in the transition; HR: Huang-Rhys factor.

Type
Transition

Carbon-on-nitrogen (CN)
The ideal h-BN structure is only modified slightly with the inclusion of a C 0 N . C 0 N in D 3h symmetry exhibits a single unoccupied defect level within the band gap, indicating that only LD transitions are possible (Fig. 3c). The 0/ − 1 CTL (ZPL for hole capture on C − N ) resides at 4.40 eV. This is higher energy than earlier calculations for bulk h-BN that obtained a value of 2.84 eV [22]. This can again be primarily attributed to the smaller mixing parameter of α = 0.31 (see Sect. III A and Sect. III D 1). It is significantly larger than measurements that located the C N acceptor at 2.3 eV above the VBM [51]. The other possible transition is via CB electron capture on C 0 N , which has a ZPL of 2.81 eV. The lattice distortion associated with these transitions is only 0.37 √ amuÅ.
The transitions on C N couple strongly to high frequency phonon modes, with the spectral function (Fig. S9) exhibiting distinct peaks at 158 meV and 185 meV. The 158 meV peak consists of a single phonon mode that has a HR factor of 0.23 (12.4% of the total HR factor of 1.88). The IPR of the 158 meV mode is 27. The 185 meV peak on the other hand consists of several modes between 182 meV and 187 meV. The largest contribution to the spectral function comes from one mode at 182 meV with a HR factor of 0. 19  spectively. The normalized emission lineshape is shown in Fig. 9a.
In order to elucidate the structural origin of the PSB, the ionic displacement due to the electronic transition can be overlaid with the phonon displacement vector for the highest frequency mode at 187 meV (Fig. 5). The coupling between lattice and electronic transition in C N is dominated by the displacement of B atoms. The B ions closest to C − N experience the largest displacement upon transition to C 0 N . However, these ions are not displaced in the phonon displacement vector so the contribution to the partial HR is essentially zero. The 6 B atoms in the next shell do not displace as much but contribute much more to the partial HR factor due to a much larger overlap with the phonon eigenvector. Finally, the B atoms in the third neighbour shell contribute the most to the partial HR factor.
The vibrational coupling can also be approximated by the 1D CC diagram in Fig. 3a. We find that the effective frequencies ω eff determined from the potential energy surface are 141 meV and 146 meV for the ground state and excited state, respectively. These frequencies translate into HR factors (Eq. 10) of 1.89 for the ground state and 1.95 for the excited state (1.88 with the generating function method). The effective frequencies are determined by the coefficient a in the fitted polynomial aQ 2 + bQ 3 , from which we find a/b > 1000 suggesting that the harmonic approximation is sound.

Carbon-on-boron (CB)
The possible charge transitions on C B have ZPLs of 4.20 eV and 3.01 eV, and the structural distortion associated with these transitions is 0.34 √ amuÅ. While the phonon spectrum for the C 0 B defect contains imaginary modes, we were unable to find lower energy structure by eigenmode following. We attribute this finding to the very small energy difference associated with a displacement of the C B atom perpendicular to the h-BN plane (see Fig. S1 of the SI). The HR factor is 1.80 and there is a well defined peak in the PSBs at 185 meV, with additional phonon replicas at higher energies. We note that the ZPLs and the spectral distribution function for C B are very similar to the ones in C N (Fig. S8, Fig. 6), which would make it very difficult to distinguish between C B and C N in a photoluminescence experiment.

Carbon-on-boron-nitrogen vacancy complex (CB-VN)
The C B -V N defect is found to have a singlet ground state in agreement with previous studies [49,52] but as other studies have already pointed out a triplet electron configuration is also possible [21,53]. For the C B -V N we considered three different structures: (i ) the triplet planar configuration (T), (ii ) the singlet planar configuration (S-planar), and (iii ) the singlet structure in which the C atom is displaced out-of-plane (S).
The singlet planar (S-planar) structure is dynamically unstable confirming previous reports [54]. The triplet structure (T), on the other hand, is dynamically stable but about 0.2 eV higher in energy than the S-planar configuration with PBE (compare Fig. S4). Finally, the singlet out-of-plane structure, which is 0.53 eV lower in energy than the S-planar configuration, is both thermodynamically and dynamically stable and thus should be the equilibrium configuration of the C B -V N defect. (A careful comparison of singlet and triplet configurations, including an assessment of the role of the exchangecorrelation functional can be found in the SI.) While the singlet state is by far the most stable configuration, it can be ruled out as a SPE source due to the large value of ∆Q =2.91 √ amuÅ. The LL transition in the triplet state, on the other hand, exhibits a much smaller lattice distortion. The ZPL, located at 1.33 eV, has a large intensity and there are multiple peaks in the PSB, occurring at energies between 28 meV and 157 meV (Fig. S8). The high-energy peak is barely distinguishable from the spectral distribution function. The HR factor for this transition is 1.72, in good agreement with previous studies on the lineshape of C B -V N in the triplet state [21]. The spectral distribution function for charged transitions on C B -V N (T) exhibit a pronounced intensity on the ZPL (see Fig. S9). However, the PSB is wide and does not exhibit a shape that can be associated with experimental lineshapes. The HR factor for these transitions are 3.2 and 3.3.

Carbon-carbon dimer
The C B -C N dimer consists of C N and C B defects located on neighboring lattice sites (Fig. 2a) and has been suggested to form at high C-doping levels [51]. The C-C bond is significantly shorter than the corresponding pristine B-N distance, and a mode with higher frequency than any pristine h-BN mode is present in the vibrational spectrum of the C B -C N dimer at 195 meV.
The neutral charge state of the C B -C N dimer is thermodynamically stable for Fermi levels above 1.17 eV making charged transitions in either the 4.1 eV or 1.6 eV-2.3 eV region possible (Fig. 2b). There are four in-gap single particle levels, two of which are occupied and two unoccupied, making a LL transition possible (Fig. 3d). The emission energy for the charged LL transition on the C B -C N dimer is 3.34 eV and dipole allowed (Fig. 7b) while the lattice distortion between the ground and excited state equilibrium configurations is 0.28 √ amuÅ. The coupling of the charge neutral emission to the vibrational degrees of freedom has been analyzed with both the 1D CC diagram (Fig. 3e) and the generating function method. From the 1D CC we find that the effective frequencies are 110 meV and 127 meV for the ground and excited state with HR factors based on the effective frequencies of 1.06 and 1.23, respectively. In comparison to the 1D CC for C N , however, the third-order coefficient carries a much larger relative weight. The ratio between the second and third-order coefficients is below 300 in both cases suggesting that anharmonic effects are not completely negligible.
The spectral function S(ω) features a significant peak corresponding to coupling to the 195 meV mode (Fig. S9), which correlates with a partial HR factor of 0.44 to be compared with the total HR factor of 1.10 as computed via the spectral function. The 195 meV mode is localized, as indicated by the small value of the IPR of 4.5, and corresponds to a stretching of the C-C bond. The computed HR factors of 1.10, 1.06 or 1.23, depending on method of calculation, agree well with the measured HR factor for the 4.1 eV luminescence of 1.3 [17].
The computed normalized lineshape compares very well with experimentally measured lineshapes [18] for the 4.1 eV emission (Fig. 9b), including both the positions of the features in the PSBs and the relative intensities between the first and second peaks. We note that the results shown are for natural carbon i.e. 12 C. For 13 C the frequency of the dominating mode at 195 meV is reduced to 191 meV.

Dissociation of carbon-carbon dimer (CN-CB)
Next, we consider the effect of spatial separation on the C N -C B defect, while the C-C distance is varied between 2.90Å and 6.32Å. An inspection of the Kohn-Sham (KS) levels shows that both LD and LL transitions are possible on the dissociated C N -C B pair (Fig. 7a). The ZPL for the LD +1/0 transition ranges from 2.0 eV to 2.7 eV (see Fig. S3; all values given here include the band edge shift between PBE and HSE06 (α = 0.6)), while the LL transitions exhibit ZPLs between 2.39 eV at a separation of 2.90Å and 1.68 eV at a C-C separation of 5.79Å, which is the longest C-C distance for which there is a well defined peak in the imaginary part of the dielectric function (Fig. 7b).
The vibrational properties of the dissociated C B -C N structures are computationally demanding to obtain due to the low symmetry of the systems. Instead, coupling to bulk phonon modes is considered by utilizing the phonon eigenvectors for the ideal system. Since the lattice distortion induced by C B and C N is small, the bulk phonon modes are expected to provide a good estimate of the emission lineshape. The spectral distribution functions S(ω) for C N computed with ideal modes and defective modes are very similar (see Fig. S2). This is likely due to the limited structural relaxation relative to the ideal structure that occurs when incorporating C impurities in h-BN. For defects such as C B -V N (T) (also shown in Fig. S2) that significantly distort the lattice, the spectral distribution computed with bulk modes is not a good approximation of the spectral distribution obtained using the modes of the defect structure. This is especially true in cases where the contribution to the HR factor mainly originates from coupling to local (small IPR), defect induced modes.
The spectral distribution functions for the LL transition on dissociated C B -C N pairs obtained in this fashion are shown in Fig. 7c for distances between 2.91Å and 5.79Å. C N -C B defect pairs couple to high frequency bulk modes for all distances considered and have pronounced peaks separated by a frequency of around 187 meV with some minor variations between the different structures.

IV. DISCUSSION
The only intrinsic defect that emerges from our analysis of the vibrational spectra as a potential SPE (V B ) has a large formation energy. Large formation energies   are also obtained for many of the best candidates that are based on extrinsic defects.
The large formation energies are consistent with the observation that pristine h-BN usually exhibits a very small concentration of color centers. To create emitters  (Table I). N B -V N can be ruled out based on the magnitude of the relaxation between the states involved ∆Q, which leads to an insensibly large HR factor [2,7]. B N can be excluded on the same premises. C B -V N (T), which has been proposed as a candidate [21], is an unlikely source due to the instability of the triplet state. The more stable singlet state, on the other hand, is unlikely to exhibit a structured emission line based on the same argument as for the N B -V N defect.
V − B has been proposed to be an optically active defect and shown to be able to host transitions in the 2 eV region before [50]. The emission lineshape of the HOMO/LUMO transition on V − B exhibits a rapidly decaying spectral weight away from the ZPL similar to many measured emission spectra in the 2 eV band (Fig. 8). Very recently scanning transmission electron microscopy (STEM) images on vacancies and multivacancies in h-BN have become available [55]. V B was identified and associated with an emission energy of 1.98 eV, which is in very good agreement with ZPLs of 1.6 eV to 1.7 eV obtained from semi-local DFT calculations. (We note that semi-local DFT calculations have been found to underestimate the LL transition energy by around 0.3 eV compared with HSE06 (α = 0.25) in the case of the NV −1 center in diamond [29] and the C B -V N (T) defect in monolayer h-BN [21].) The computed lineshape for the LL transition on V − B is compared in Fig. 8a with the lineshape from Ref. 55 that has been associated with V B . While some features agree between computed and experimental lineshapes such as the position of the first peak close to the ZPL and the peak at around 162 meV, which is more pronounced in the case of emission from D 3h , the overall agreement is poor. Specifically, the ZPL intensity is much larger in the computation and the region between the first peak and second peak carries a large spectral weight in the computation while there is a pronounced gap in most measurements.
While one could question whether semi-local DFT calculations are sufficient for modeling excited state geometries of V − B , the comparison with the 2.25 eV emitter measured in Ref. 7 (Fig. 8b) suggests that the lineshape observed in Ref. 55 rather originates from some other defect. Specifically, the computed lineshape for the V − B in D 3h geometry is in good agreement with the measured spectrum, which supports the assignment of in-gap transitions on V − B to at least some emitters in the 2 eV region. The ZPLs for isolated carbon impurities (C N , C B ) as well as the nearest-neighbor C B -C N dimer fall outside of the 1.6 to 2.3 eV window considered here. Nextnearest and farther neighbor C B -C N pairs (C-C distance ≥ 2.90Å) exhibit, however, varying ZPLs around 2 eV depending on separation distance (Table II). Unlike the (nearest-neighbor) C B -C N dimer, there are no direct C-C bonds present in these configurations and the PSBs mainly originate from coupling to bulk modes. The lineshape exhibits only small changes with increasing C-C distance (Fig. 7) while the ZPL varies strongly (Table II), which is a key feature of the SPEs found experimentally in the 2 eV region.
The present analysis demonstrates that dissociated C B -C N defects correspond to a range of different ZPLs with very similar PSBs. In practice, stabilization of these different ZPLs requires confinement of the atoms at specific lattice sites. Since h-BN is a strongly covalent material bond breaking is energetically costly and migration barriers are high. As a result, it is plausible that C B -C N defects can be stabilized at a range of distances. The computed lineshapes for C N and C B show excellent agreement with the lineshape of a 4.1 eV emitter reported in Ref. 17 (Fig. 9a). The frequencies of the dominant modes in C N and C B are 182 meV and 187 meV, in excellent agreement with the measured frequencies. It is, however, important to note that neither defect distorts the lattice significantly enough to induce local modes, which was argued in Ref. 17 to be the origin of the PSB. We note that the excellent agreement might be partly coincidental since in bulk samples, such as the one measured in Ref. 17, the longitudinal optical-transverse optical (LO-TO) splitting is significant and might push the PSB away from the ZPL by ∼ 10 meV. In fact, the highest intensity peak in the PSB for the 4.1 eV emitter has been suggested to originate from coupling to the zone center LO mode at 200 meV [20]. The present calculations for C N and C B yield lineshapes and ZPLs consistent with the 4.1 eV emitter observed experimentally.
Our calculations also suggest C B and C N to be strong candidates for emitters found at 3.2 and 3.4 eV [19,57]. The lineshapes of both C N and C B compare well with the 3.22 eV emitter in Ref. 19, where the first PSB was found at 200 meV from the ZPL with two additional distinct phonon replicas (Fig. 6). The experimental data was recorded at room temperature, which explains the considerably broader spectrum compared to the calculation. The position of the first PSB and the phonon replicas agree well between experiment and calculation safe for a ∼ 10 meV offset, which, as noted above, arises from the presence of LO-TO splitting in a bulk sample. Both ZPLs and lineshape thus suggest that C N and C B can be associated with an 3.2 eV emitter.

Carbon-carbon dimer (CB-CN)
As noted above, the 4.1 eV emission has also been associated with the (nearest-neighbor) carbon-carbon dimer (C B -C N ) [56]. While DFT calculations based on semilocal exchange correlation (XC) functionals yield a ZPL of 3.34 eV for this defect, hybrid functionals with a larger fraction of exact exchange predict a ZPL in much better agreement with the measured 4.1 eV emission line [56]. Here, it is important to keep in mind that the LL transition involving C B -C N should be much less affected by the band gap error from the LD transitions on isolated C N or C B since it only involves localized states, which are already well described by semi-local XC functionals.
Our analysis shows that the lineshape of the C B -C N defect agrees very well with a measured emission spectrum from Ref. 18 (Fig. 9b). The PSBs are located at approximately the correct positions, although the spectral weight of the ZPL is slightly larger for the computed lineshape. Interestingly, the 195 meV mode that is the origin of the PSB in the C B -C N dimer is indeed a local mode in contrast to the dominant modes in C N and C B , which are bulk-like. Furthermore, the computed HR factor of 1.1 to 1.2 compares favorably with a measured value of 1.3 for the 4.1 eV emitter [17].
The Stokes shift on the excited state landscape of 0.2 eV obtained with HSE06 (α = 0.25) is notably larger than the value of 0.13 eV from PBE. A larger Stokes shift would indicate a smaller spectral weight of the ZPL. Assuming that this difference is dominated by the displacement of the potential energy surfaces relative to each other, the spectral function can be renormalized and, in fact, renormalizing the electron-phonon spectral function to a HR factor of 1.54 leads to an even better agreement.
We note that since the PSB originates to a large extent from the local mode at 195 meV, which predominantly involves C motion, isotopic effects may appear in the position of the PSB. This is in contrast to, e.g., the transitions on C N and C B , where isotopic control over the C atoms in the host matrix should not affect the PSB. Hence if one can isotopically control the formation of the C B -C N dimer to comprise a pair of 13 C atoms instead of the naturally occurring 12 C one should detect a small variation in the position of the PSB if the 4.1 eV line originates from the C B -C N dimer. Based on our calculations we estimate this frequency shift to be 4 meV, placing the local C-C mode at 191 meV for 13 C. While emitters based on 13 C have been fabricated [58], changes in the PSB have not been investigated yet.
The transitions on intrinsic defects that have ZPLs in the vicinity of 4.1 eV include hole capture on V 0 N , as well as hole capture on V − B , which give rise to ZPLs of 4.04 eV and 4.12 eV, respectively (see Table I). However, the structural distortions associated with these charged (LD) transitions are relatively large and unlikely to be compatible with the measured lineshape and HR factor (Table I).

V. CONCLUSIONS
We have analyzed the optical transitions occurring from charge transitions (LD) and charge neutral (LL) transitions for a wide range of defects in h-BN. For a selection of these defects, we have examined the possibility of assigning defect emission to well known SPEs by calculating line shapes and emission spectra with the generating function method. We have found that V B , C N , C B , and C B -C N can host transitions that couple strongly to high-frequency modes while still exhibiting a moderate HR factor between 0.9 and 2.5.
The main conclusions are (i ) V B −1 is likely to be a SPE with a narrow emission band with a ZPL in the 1.6 eV to 2.3 eV region. The lineshape shows decent agreement with the measured lineshape of an emitter with a 2.25 eV ZPL.
(ii ) The lineshapes of C N and C B are in excellent agreement with measured lineshapes for the 4.1 eV emitter. Furthermore, the ZPLs are in good agreement with the 4.1 eV ZPL. In addition, we note that these defects also exhibit ZPLs for the reverse transition that agrees well with the observed 3.2 to 3.4 eV luminescence.
(iii ) Next-nearest and farther neighbor (dissociated) C B -C N defect pairs allow for both charge neutral and charged transitions over a wide range of ZPL energies with narrow emission bands.
(iv ) Our findings support the assignment of the 4.1 eV emission to the C B -C N dimer (nearest-neighbor pair) based on the lineshape but we note that there are additional defects that exhibit similar emission lineshapes.
The current results corroborate the emerging consensus that the low frequency band (1.6 to 2.3 eV) originates from V B based defects and the high frequency (∼ 4.1 eV) emitter originates from C based defects. We note that recently additional defect configurations have been suggested that have not been considered here, including Stone-Wales defects [59] as well as oxygen-carbon pairs [60].