Phonon-assisted emission and absorption of individual color centers in hexagonal boron nitride

Defect centers in hexagonal boron nitride represent room-temperature single-photon sources in a layered van der Waals material. These light emitters appear with a wide range of transition energies ranging over the entire visible spectrum, which renders the identification of the underlying atomic structure challenging. In addition to their eminent properties as quantum light emitters, the coupling to phonons is remarkable. Their photoluminescence exhibits significant side band emission well separated from the zero phonon line (ZPL) and an asymmetric broadening of the ZPL itself. In this combined theoretical and experimental study we show that the phonon side bands can be well described in terms of the coupling to bulk longitudinal optical (LO) phonons. To describe the ZPL asymmetry we show that in addition to the coupling to longitudinal acoustic (LA) phonons also the coupling to local mode oscillations of the defect center with respect to the entire host crystal has to be considered. By studying the influence of the emitter's wave function dimensions on the phonon side bands we find reasonable values for size of the wave function and the deformation potentials. We perform photoluminescence excitation measurements to demonstrate that the excitation of the emitters is most efficient by LO-phonon assisted absorption.

Abstract. Defect centers in hexagonal boron nitride represent room-temperature single-photon sources in a layered van der Waals material. These light emitters appear with a wide range of transition energies ranging over the entire visible spectrum, which renders the identification of the underlying atomic structure challenging. In addition to their eminent properties as quantum light emitters, the coupling to phonons is remarkable. Their photoluminescence exhibits significant side band emission well separated from the zero phonon line (ZPL) and an asymmetric broadening of the ZPL itself. In this combined theoretical and experimental study we show that the phonon side bands can be well described in terms of the coupling to bulk longitudinal optical (LO) phonons. To describe the ZPL asymmetry we show that in addition to the coupling to longitudinal acoustic (LA) phonons also the coupling to local mode oscillations of the defect center with respect to the entire host crystal has to be considered. By studying the influence of the emitter's wave function dimensions on the phonon side bands we find reasonable values for size of the wave function and the deformation potentials. We perform photoluminescence excitation measurements to demonstrate that the excitation of the emitters is most efficient by LO-phonon assisted absorption.

Introduction
Single-photon emitters are at the heart of many promising quantum technologies such as quantum computing and quantum cryptography. Although there are various solid-state emitters of single photons, there is still no system available, which simultaneously meets all requirements [1]. For example, self-assembled semiconductor quantum dots [2,3] represent a mature technology platform, but presently yield optimum performance only at cryogenic temperatures. Another class of single-photon emitters are defect centers in insulators, such as color centers in diamond [4,5,6]. These defect centers exhibit prominent single-photon emission also at room temperature but suffer from a high refractive index and challenging processability of the host crystal. Recently, a new class of single-photon emitters in atomically thin semiconductors has been discovered [7,8,9,10,11,12,13] and has been deterministically positioned on the nanoscale by strain engineering [14,15,16]. Defect centers in hexagonal boron-nitride (hBN) combine features from these classes [17,18,19,20,21,22,23,24,25,26,27,28,8,30,31,32]. They have the characteristics of an atomically sized defect center and at the same time share the advantages of a layered structure, i.e., the ultimate limit of miniaturization due to their atomic thickness and high mechanical robustness.
In this work, we investigate the photoluminescence (PL) and photon absorption of defect centers in hBN nanocrystals. By developing a theoretical model to calculate the PL spectrum taking into account the coupling to bulk longitudinal optical (LO) and longitudinal acoustic (LA) phonons and to the oscillation of a local mode, we explain measured emission spectra achieving an excellent agreement between experiment and theory. In addition, we verify by means of photoluminescence excitation (PLE) spectroscopy that coupling to LO phonons provides an efficient way of exciting individual single-photon emitters. This then allows for the isolation of a desired single-photon emitter from the wide range of emission energies of various defect centers in hBN.

Emission spectra
Localized single-photon emitters are investigated in hBN nanopowder deposited on a Si/SiO 2 substrate (see Supplementary Material for details about the structure). Figure 1 presents typical measured room-temperature PL spectra of six different localized light emitters in hBN. The transition energy of the emitters E ZPL , which is called the zero phonon line (ZPL) (labeled with numbers (1) to (6) in Fig. 1(a)) varies over a large range between 1.6 eV and 2.5 eV [19]. The energy of the excited states of the emitters has already been extensively studied, both theoretically via DFT calculations [33,34,35,31] and experimentally [19,23,24,25,26,27,8,31,32,36]. Although different types of atomic defects such as N-or B-vacancies or substitutions with carbon or oxygen might exist, the distribution of transition energies is remarkably homogeneous. Therefore, it is unlikely that each emission energy in hBN stems from a different type of defect. One promising proposal for the measured wide range of transition energies is that at least part of the energy spread is due to a Stark shift [31] resulting from an electric field which arises from a trapped charge near the color center [23]. Another suggested origin is a local strain distribution in the crystal [37,28,30,36]. Our results suggest that both effects might contribute to the broad energy spread of the ZPL energies.
The narrow ZPLs of the emitters are accompanied by characteristic phonon side bands at lower energies. To compare the phonon side bands of the different emitters more easily, the spectra are plotted against the detuning from the ZPL energy in Fig. 1(b). Two side bands well separated from the ZPL between 0.2 eV and 0.15 eV below E ZPL are visible, which can be attributed to two longitudinal optical (LO) phonons [21]. Also the ZPL itself is always more or less asymmetric with a steeper high energy side, which hints towards the contribution of low-energy (LE) phonon modes with energies below 50 meV. As will be shown later, we also assume two types of LE modes to be important. Comparing the phonon side bands of the different emitters in Fig. 1(b) we find significant variations. Not only the relative height of the LO phonon side bands compared to the ZPL intensity, but also the detailed structure, i.e., the relative heights of the LO 1 and LO 2 peak with respect to each other, vary from emitter to emitter. Having a closer look at the LE phonon modes near the ZPL, also this spectral phonon feature has different shapes and strengths for different emitters as highlighted in Fig. 1(b). In the following sections we will study these phonon side bands and develop a detailed understanding of their origin. In this context we will discuss which arguments support the identification of the LE phonons as longitudinal acoustic (LA) phonons and local mode oscillations of the defect with respect to the entire crystal.

Phenomenological analysis
The properties of the LE and LO phonon side bands change significantly among the emitters. Therefore, we first try gaining a broad, quantitative overview of the different features. For this purpose we use a phenomenological model to fit the room temperature PL spectra of 165 different emitters. Figure 1(c) shows an exemplary spectrum (blue) and its phenomenological fit (orange). It is obvious that two LO phonon modes with energies at E LO,1 ≈ 165 meV and E LO,2 ≈ 195 meV have to be taken into account (red line). These energies are also supported by Raman measurements and ab initio calculations [1]. While trying different models for the LE phonons we found that at least two LE modes are required to achieve reasonable fits for all measured spectra. The best agreement was found for E LE 1 = 14 meV and E LE 2 = 30 meV. We will later show that these energies can be well identified as representative for LA phonons in the case of LE 1 and for a local mode in the case of LE 2 . But for now we will stick to the label LE phonons for the low energy side bands. The applied fit function for the PL intensity reads (3)   (6). (b) Same as (a) but plotted as a function of the detuning with respect to the ZPL energy to highlight the different shapes of the phonon side bands. (c) Phenomenological fit (orange) of a measured PL spectrum (blue). The constituents of the fit are the ZPL (green), single phonon LO side bands (red), the low energy (LE) side bands LE 1 (yellow) and LE 2 (violet) with the respective absorption peaks (dotted) and the two phonon side bands (dashed turquoise).
where we include a ZPL, single-phonon processes for the LO phonons and single-and two-phonon processes for the LE phonons. We also take single-phonon absorption processes for the LE modes into account. For all peaks we assume a Lorentzian shape.
From each fit we retrieve the integrated peak weights A j (j = ZPL, LO 1 , LO 2 , LE 1 , LE 2 ) and respective peak widths γ k (k = ZPL, LE 1 , LE 2 ), where we assume that the LO side bands have the same broadening as the ZPL, which is a good approximation for dispersionless LO modes. Details about the fit function can be found in the Supplementary Material. In Fig. 1(c) we additionally show the two LE side bands (yellow for LE 1 and violet for LE 2 ) in phonon emission (solid) and absorption (dotted). The dashed turquoise line depicts the contribution from two phonon processes. Note that the absorption and two-phonon side bands are strictly related to the corresponding singlephonon emission side bands and therefore do not introduce additional fit parameters. The fit slightly deviates from the measurement for energies larger than the ZPL energy. The reason is, that a longpass filter is used close to the transition energy to block the exciting laser. This often reduces the measured high energy side of the ZPL. Therefore, we only used energies smaller than E ZPL for the fits.
For the large data set of 165 emitters we identify correlations between the different spectral features. We start by comparing the LO phonon side bands to the ZPL. For this purpose, Fig. 2(a) shows the LO contribution, i.e., the ratio between the LO side band and the ZPL weight, as a function of the ZPL energy. The colors of the dots represent the energy of the driving laser. The photon energies of the three lasers are marked by the dashed lines. We find transition energies almost equally distributed between 1.6 eV and 2.2 eV [22,23,25,28,36]. Even a few emitters appear above 2.4 eV. Larger transition energies could not be measured, because of the longpass filter blocking the exciting laser. Having a look at the values of the LO contribution, i.e., the vertical distribution of the dots, we find that they have a clear tendency following a positive slope. While for transition energies below 1.9 eV most of the points lie between 0 and 0.2, they appear between 0.1 and 0.45 for energies around 2.1 eV. This shows that the emitters with smaller transition energies couple on average less efficiently to the optical phonons. This is in good agreement with the assumption that the shift of the ZPL energy is associated with the Stark effect. A static electric field acting on the emitter's dipole should increase the distance between the different charge centers. In consequence this should also influence the coupling to the polar LO phonons. The additional spread of LO side band strengths for a fixed energy, especially around 2.1 eV, could be due to the effect of local strain. Distortions of the lattice constants in the vicinity of the defect center should not only influence the distance between the charge centers, but also change the size of the emitter's wave function and therefore the effective distance between the charge centers. Later on we will discuss the influence of these properties on the phonon coupling in more detail.
To shed more light on the properties of the LO phonons, Fig. 2(b) shows the energies of the LO phonon side bands as a function of the ZPL energy. We find two clusters of data points, which we associate to the two LO phonon modes near the Γ point. The energies of LO 1 (yellow) appear around 165 meV, those of LO 2 (violet) around 195 meV, which agrees well with calculated band structures and Raman measurements [1]. As indicated by the dotted lines, we find that the phonon energies have a slight trend to   higher values for larger ZPL energies. As we identify these energies with the phonon energies of the hBN bulk crystal, this correlation suggests that the origin of the energy shift of the ZPL should also influence the LO phonon energies. This could be the case for a local strain field, which naturally changes the phonon dispersion. But also the dielectric environment might alter the LO phonon energies, as it affects in particular the LO-TO splitting [39]. Figure 2(c) focusses on the LE phonon side bands by plotting the ratio between the first side band weight A LE 1 and the ZPL weight A ZPL as a function of the ZPL energy. The dots are scattered over a wide range from 0 to almost 1 for any transition energy. However, the majority of the points is equally distributed between 0.1 and 0.7. This finding does not provide any strong correlation between the coupling to the LE 1 (LA) phonons and the transition energy. The color and size of the dots represents the LE 1 contribution to the entire LE side band, i.e., A LE 1 /A LE , where A LE = A LE 1 + A LE 2 (small dots are blue, large dots red). Also this quantity shows no correlation in the dot pattern. Small blue and large red dots are found everywhere. This suggests that the coupling to LE phonons is not strongly affected by the parameters that are expected to govern the transition energy, i.e., the distance of the charge centers. However, we find a wide spread of dots on the vertical axis, i.e., a non-negligible variation of coupling strengths for a given transition energy. This finding will be traced back to variations of the emitter's size as discussed in more detail later.
Finally, Fig. 2(d) presents the ZPL width as a function of ZPL energy. We find a slight trend to smaller line widths at small transition energies, which is a similar to what has been found in Fig. 2(a) for the strength of the LO side bands. This is in line with established models which associate the ZPL width with second-order phonon coupling mechanisms [40]. The mechanism suggested in Ref. [40] relies on the scattering of acoustic phonons with the polaron formed by optical phonons. Consequently, a weaker coupling to LO phonons should result in a smaller lattice displacement of the polaron and therefore in a smaller scattering rate for the acoustic phonons. This could lead to the narrowing of the ZPL we found here. Another reason for this trend could be a stronger spectral wandering induced by trapped charges which are also considered as the origin of the Stark shift of the ZPL energies [23].

Theoretical model
The phenomenological model used to fit the 165 spectra in Sec. 2.2 already includes some assumptions about the expected phonon modes, but it does not allow to retrieve further information about the emitters themselves. To support our assumptions about the mode energies we use a microscopically motivated model to reproduce the measured PL spectra. We also use this model to find possible origins for the correlations and non-correlations of the different phonon features and the ZPL energy found in Fig. 2. As already discussed in detail, our model should include the coupling to LO phonons to reproduce the two side bands around 165 meV and 195 meV below the ZPL. Also the coupling to acoustic phonons, which typically leads to a low-energy broadening of the ZPL, should be present in these systems. From defect centers in diamond, such as the extensively investigated NV − center, it is well known that the prominent phonon side band is dominated by the coupling to local mode oscillations [41,42]. These modes describe displacements of the defect atoms with respect to the entire host crystal. Due to the similarity of the atomic defect structures, we expect that also in hBN local modes might have a significant influence on the PL spectrum.
We model the PL intensity spectrum of an emitter by calculating the optical susceptibility for a two-level system coupled to phonons, which in the time domain is given by [43,44,45,46] where θ(t) denotes the Heaviside step function, ω 0 = E ZPL is the polaron-shifted transition energy of the emitter, i.e., the energy of the zero phonon line, T 2 the dephasing time and Φ describes the phonon-induced dephasing with This phonon influence is determined by the temperature T and the phonon spectral density J j at the phonon energy ω ph for the phonon mode j. Note that such a model was commonly used to simulate the phonon coupling of F-centers [43]. In this section T = 300 K is considered to allow comparison with the experimental results in Sec. 2.2.
In the case of coupling to bulk phonons the spectral density is given by [47] J Here, g j (q) denotes the coupling matrix element for the coupling to a phonon with wave vector q in the phonon branch j and ω j (q) is the corresponding phonon dispersion relation. With Eq. (1) the absorption spectrum α(ω) of the emitter is simply given by the imaginary part of the Fourier transform From this expression the emission spectrum I(ω) is then retrieved by inverting the absorption α with respect to the ZPL at ω 0 , i.e., We account for the coupling of the defect to LO 1 and LO 2 phonons via the Fröhlichinteraction and for the coupling to LA phonons via the deformation potential coupling where V is a normalization volume, is the mass density, ε s and ε ∞ are the static and high frequency dielectric functions, respectively, and D e and D h are the deformation potentials for electrons and holes, respectively. The form factors F e,h (q) are determined by the wave function of the excited state of the emitter as we explain in the following. Another common coupling mechanism between charges and LA phonons is via the piezoelectric effect. In Ref. [48] it was shown that the piezoelectric constants vanish with an increasing number of hBN layers. As we are investigating nanopowder samples, we are dealing with small bulk crystals with large layer numbers. The thickness of the investigated hBN crystals is continuously distributed between 40 nm and 100 nm as exemplarily shown in the Supplementary Material. In Ref. [48] it was demonstrated that due to inversion symmetry the piezoelectric constants vanish for even layer numbers. If this coupling mechanism to LA phonons played an important role, it should only appear for half of the studied defects, i.e., those located in samples with odd layer numbers. This effect should lead to a separation of data points for crystals of even and odd layer number in Fig. 2(a) but one does not find such a feature. Therefore, we conclude that the piezoelectric coupling to LA phonons does not play an important role in our samples and assume that this coupling mechanism can be neglected.
One central aspect of hBN is its pronounced anisotropy arising from its layered structure. In our model we take this into account by distinguishing in-plane and out-ofplane directions. We assume that the dipole of the emitter lies in the plane of one hBN layer. Therefore, we consider a wave function of the defect's excited state consisting of two differently charged centers with an in-plane distance d as schematically shown in Fig. 3(a). The negative charge will be referred to as electron (e) and the positive one as hole (h). Each of the charge centers has a Gaussian wave function with an in-plane localization length a and an out-of-plane localization length b. This form of the wave function can be seen as an approximation for a wave function typically calculated from ab initio theory. The advantage of our approach is the flexibility to easily study the influence of changes of geometrical quantities. With this form of the wave function the form factors in Eqs. (6a) and (6b) read where the phonon wave vector q is split into an in-plane q r and an out-of-plane q z component and d is the vector connecting the two charge centers. After integrating (3) over the in-plane angle of the phonon wave vector, the influence of the wave function geometry on the coupling strengths is given by and where J 0 is the Bessel function of first kind and of zeroth order. Due to these form factors, coupling to phonons is essentially restricted to the wave vector region q r a −1 , For the Fröhlich coupling, different dielectric constants for the in-plane and the out-of-plane direction are considered. The energies of the phonons, which enter the model, are extracted from the simulated dispersion relation shown in Fig. 3(b) adapted from Ref. [1]; we assume two LO modes (LO 1 and LO 2 ) with constant dispersions. For the LA phonons we consider linear dispersions, where we again distinguish between the in-plane and the out-of-plane direction. More details about the model parameters can be found in the Supplementary Material. As will be shown below, the 14 meV LE 1 mode introduced in Sec. 2.2 can be seen as a reasonable approximation for the side band stemming from the LA phonons. In Ref. [34] it was calculated that the C B N V center has the strongest coupling to an in-plane breathing mode of the defect center, which has an energy of 30 meV. Also other studies suggest this atomic structure as promising candidate for the hBN color center [49,50]. This motivates us to interpret the LE 2 mode in the phenomenological model in Sec. 2.2 as representative of such a breathing mode. The coupling to this local mode is modeled by a Lorentzian-like distribution of the phonon spectral density [51,52,53] J LOC (ω ph ) = 1 π The parameters entering this coupling are the local mode energy ω LOC = E LOC , the coupling strength g and the width ∆. For the considered parameters the spectral shape of the side band due to the local mode agrees well with the ordinary Lorentzian shape considered for the phenomenological fits in Sec. 2.2. We show in the Supplementary Material that our model for the local mode coupling reproduces PL spectra of different defect centers in diamond very well, showing the accuracy of the model. To reproduce the finite width of the ZPL, we assume a constant dephasing rate 1/T 2 , providing us with a dephasing time T 2 . We note that this dephasing time includes pure dephasing, homogeneous broadening mechanisms as well as inhomogeneous processes, e.g., spectral wandering [27]. While the phonon dispersion relations and dielectric constants are well known for hBN, the deformation potentials are only hardly known. Simulations indicate that the relevant values are in the range of some hundred meV to a few eV [54]. We therefore will try finding reasonable values for the deformation potentials D e and D h for the LA phonon coupling of the considered defect centers in hBN. The other unknown parameters in our model are the in-plane and out-of-plane localization lengths a and b, respectively, and the charge separation d in the defect state. To be able to extract information on these parameters from the measured spectra, we will analyze their role for the coupling to LO and LA phonon modes, respectively.
For the coupling to LO phonons the only unknown quantities are the geometrical parameters in-plane localization length a, out-of-plane localization length b, and charge distance d. To investigate their influence, we first concentrate on the coupling to LO phonons. In Fig. 4(a) we plot the contribution of the LO single-phonon side band, i.e., the fraction between the weight of the first LO phonon side bands of LO 1 and LO 2 and the ZPL weight, as a function of the charge distance d for different localization lengths a between 0.1 nm (blue) and 1 nm (red). The distance d is varied on the same scale. This should be a reasonable range for these quantities as the in-plane lattice constant of hBN is 0.4 nm and the extension of an atomistic defect state should be of similar size. Under the assumption that the Stark effect is responsible for the different transition energies of the emitters, we additionally plot the measured data points from 2(a) as gray dots into Figs. 4(a). The relation between the energies of the experimental values of Fig. 2(a) (gray dots in Fig. 4(a)) and the charge distances d in Fig. 4(a) is obtained by the fits of the measured spectra in Sec. 2.4, which yield values between d = 0.2 nm and d = 0.5 nm. By the given choice of the ZPL energy axis the majority of the measurement points appear in this range. The simulated LO contributions vanish for vanishing in-plane distance d because in this limit there is no dipole moment associated with the defect. For a given in-plane localization length a the LO contribution grows with increasing d. When increasing the localization length a the curves rise slower, because the Fröhlich coupling strength depends on the effective dipole strength of the wave function. To achieve the same strength for larger localization lengths a the charge separation has to be accordingly larger. The increase of the LO contribution with increasing distance nicely resembles the distribution of the measured dot pattern, which  supports the interpretation that the shift of the ZPL is at least partially related to a Stark shift. To further analyze the influence of the LO phonons, we show exemplary normalized PL spectra for different values of a and d in Fig. 4(b). The colors agree with the corresponding curves in (a). At the bottom for a = 0.5 nm the two single-phonon LO side bands are well resolved around −0.2 eV, also the three two-phonon side bands between −0.3 eV and −0.4 eV can be seen. As expected the amplitude of the side bands increases for larger in-plane distance d. At the top of Fig. 4(b) the same spectra are shown for a = 0.1 nm. Here, the LO side bands grow very rapidly when increasing d, they even become larger than the ZPL. For the largest considered d the two-phonon side bands are even larger than the single-phonon side bands. Indeed, from Eqs. (2) and (3) it can be shown that, as long as successive side bands do not spectrally overlap, the contributions of single-phonon (A LO ) and two-phonon (A 2LO ) side bands with respect to the ZPL are given by Thus, as soon as A LO > 2A ZPL the two-phonon sideband will exceed the single-phonon one. However, in the measured PL spectra of hBN color centers this regime where the side bands exceed the ZPL is not reached. Therefore, we can conclude that the localization length should lie between a = 0.3 nm and 1 nm and the charge distance between d = 0 nm and 1 nm. After finding reasonable values for the in-plane size of the wave function (a and d), in Figs. 4(c) we focus on the values for the deformation potentials D e and D h . When considering the special case of a vanishing charge distance, d = 0, the deformation potential coupling strength in Eq. (8b) is proportional to (D e − D h ) 2 . So we will use this limit to find reasonable values for the difference of the deformation potentials. To provide a quantitative measure, we calculate the LA contribution to the entire spectrum via This value is plotted as a function of |D e − D h | in the inset in Fig. 4(c) for wave function dimensions of a = 0.4 nm, b = 0.2 nm and d = 0 nm. For D e = D h the LA contribution vanishes and no LA side band appears in the spectrum, while the side band dominates the entire spectrum for large values of |D e − D h |. In Fig. 2(c) we found that the LA contribution in the measured spectra also for small transition energies, i.e., small d, yields values between zero and one. Therefore, reasonable values for the difference of the deformation potentials should lead to an LA contribution of around 0.5. We find that this is fulfilled for |D e −D h | ≈ 0.6 eV. However, the deformation potentials not only determine the LA contribution but are also responsible for the shape of the spectra. To demonstrate this influence, we show exemplary spectra for different values of |D e − D h | still for a vanishing dipole moment d = 0 in Fig. 4(c). We find that with increasing |D e −D h | (blue to red) the asymmetry of the ZPL peak increases till the entire spectrum is dominated by the LA side band for the 1 eV case (red). Here the ZPL is only a rather small peak on the broad phonon background. From here on we fix |D e − D h | = 0.6 eV.
As already mentioned, we assume that the coupling to the LE 1 mode phenomenologically describes the coupling to LA phonons. One might wonder why the coupling to a discrete mode with a Lorentzian lineshape approximates reasonably well the coupling to the continuum of LA phonons. To confirm that this is indeed a valid approximation, the black solid line in Fig. 4(d) shows a simulated room-temperature PL spectrum, where the coupling to the continuum of LA phonons is fully taken into account. This spectrum was then fitted by the phenomenological model of Sec. 2.2 yielding the orange line. The contributions of the two LE modes (yellow and violet) and the ZPL (green) in the phenomenological model are plotted separately as before. We find an excellent agreement between the simulated spectrum including only LA phonons and the phenomenological fit. Furthermore we find that only the LE 1 mode has a relevant contribution in the fit. This shows that indeed the LE 1 mode is a reasonable approximation for the contribution of the LA phonons to the spectrum.
In the next step we investigate the influence of the sum of the deformation potentials. For a non-vanishing charge distance d the Bessel function in Eq. (8b) does not vanish leading to the fact that the phonon coupling depends on D e and D h separately. In Fig. 5(a) we show simulated PL spectra for a = 0.4 nm, b = 0.2 nm and d = 0.5 nm and different values of D e +D h between 0 eV (blue) and 1 eV (red). We find that the weight of the LA side bands close to the ZPL strongly increases with the sum of the deformation potentials. It is also nicely visible that the LO phonon side bands are additionally broadened by the LA phonons, demonstrating that multi-phonon contributions are well reproduced by the model. In Fig. 2(c) we found that the LA contribution (LE 1 mode) to the spectrum exhibits no significant correlation with the transition energy, i.e., with d. This finding should be reproduced by the choice of D e + D h . Therefore, in Fig. 5(b) we plot the LA contribution as a function of the charge distance d for different values of D e + D h between 0 eV (blue) and 4 eV(red) (a = 0.4 nm and b = 0.2 nm). All curves start at the same value of slightly below 0.6 for d = 0 nm. We find that for deformation potential sums of 1 eV or larger the LA contribution grows significantly with increasing d. One would therefore expect a clear correlation between these two quantities. This makes us conclude that D e + D h should be smaller than 1 eV. In the special case of D e + D h = 0.6 eV = |D e − D h | the LA contribution is a straight line, because in this case either D h = 0 eV or D e = 0 eV, which makes the charge distance irrelevant. However, we do not expect one of the deformation potentials, either for electrons or for holes, to vanish in hBN. This would mean that either electrons or holes do not scatter with LA phonons via the deformation potential coupling, which makes us disregard this value. We found good results for D e + D h = 0.2 eV (D e = 0.4 eV, D h = −0.2 eV) as will be demonstrated below in Sec. 2.4 when directly comparing the simulations with measured PL spectra.
After having selected reasonable values for the deformation potentials, we now study the influence of the other two geometrical parameters, the in-plane localization length a and the out-of-plane localization length b of the wave function. Figures 5(c,d) focus on the in-plane localization length a. In Fig. 5(c) exemplary spectra for d = 0.5 nm, b = 0.2 nm and different values of a between 0.2 nm (blue) and 1 nm (red) are shown. While the LO phonon side bands increase significantly when reducing a, as discussed earlier, the LA phonon side band next to the ZPL does not change significantly. To provide a quantitative measure for the influence of a on the LA coupling in Fig. 5(d) we plot the LA contribution as a function of a for different values of d. All curves slightly decrease from values between 0.6 and 0.7 to 0.4 for a = 1 nm. This spread can give rise to part of the spread of LA contributions found in Fig. 2(c). Finally, the influence of the out-of-plane localization length b is considered. In Fig. 5(e) we again show PL spectra for a = 0.4 nm, d = 0.5 nm and different values of b between 0.2 nm (blue) and 1 nm (red). Two pronounced features are found in the spectra: (i) the LA phonon side band increases with decreasing b; (ii) the two LO single-phonon side bands LO 1 and LO 2 slightly increase with decreasing b, but not by the same rate. While LO 1 and LO 2 are almost equally strong for b = 1 nm (red), LO 1 is significantly stronger for 0.2 nm (blue). These two features are quantified in Fig. 5(f). The solid lines show the LA contribution and the dashed lines the ratio between the peak intensities of the two LO single-phonon side bands I LO 2 /I LO 1 as functions of b for different values of a. The LA contribution rapidly drops from values up to 0.9 for b = 0 nm to 0.2 for b = 1 nm. This spread of LA phonon side band contributions agrees well with the measured values found in Fig. 2(c). The fraction of the LO side band peaks increases for growing b. This variation of the shape of the LO side band was also found in the measured spectra in Fig. 1. The reason for this variation of the LO side bands is that the larger LO phonon energy is not present for the out-of-plane phonon wave vectors, as can be seen in Fig. 3(b). So b mainly influences the LO 1 phonon side band. A larger out-of-plane localization length b in Eq. (8a) leads to a smaller range of out-of-plane phonon wave vectors that is excited by the emitter, which leads to a weaker effective coupling strength. Therefore the fraction I LO 2 /I LO 1 reduces for growing b. The out-of-plane lattice constant of hBN is 0.66 nm, which makes the considered values for b a reasonable range for this parameter.
This parameter study shows that slight variations of the geometrical parameters of the defect's wave function in the range of reasonable dimensions lead to strong changes in the PL spectra. In the next step, we apply this model to directly reproduce measured PL spectra by varying the geometrical parameters. We also take the local mode for the low-energy coupling in the range of 30 meV into account and vary its coupling parameters.  Figure 6(e) depicts a temperature series taken on the same emitter. The temperature decreases from top to bottom from 300 K to 20 K. In addition to the experiments (colored lines) we show corresponding simulations as black lines. All parameters for the simulations are listed in Tab. 1.

Comparison between experiment and theory
Comparing the different spectra, we find a convincing agreement between the measured and the calculated spectra. The theoretical approach used in Ref. [24], which considers the phonon spectral function as a free parameter, also led to good agreements with measured spectra. However, the authors did not draw conclusions on the emitters' wave functions and details of the electron-phonon coupling matrix elements. From the parameters for the dimensions of the wave functions in our model listed in Tab. 1 we find that all sizes (a, b, and d) are in the range between 0.1 nm and 0.5 nm. This is a reasonable parameter range, because it shows that the assumed emitter wave function spreads over a few unit cells at most and is mainly restricted to the layer in which the defect is located. While the LO phonon energies scatter slightly around the already  discussed values of E LO 1 ≈ 165 meV and E LO 2 ≈ 195 meV, the energy of the local mode is rather fixed. This might have two reasons. On the one hand, the side band in the PL spectrum partly overlaps with the LA contribution at room temperature, which makes a precise determination of the peak position less easy. On the other hand, strain distributions could have less influence on the energetic position of the mode, but rather change the coupling strength g and the width ∆. We want to remark that, the low temperature spectrum in Fig. 6(d) shows a pronounced peak at E ≈ 2.05 eV, i.e., 35 meV below the ZPL. The dashed line includes this peak as a local mode in the simulation. However, this spectral feature could also stem from another nearby emitter. Usually, PL spectra at cryogenic temperatures exhibit a large number of narrow lines in these samples, which makes a clear identification of the lines difficult. Table 1. Model parameters for the simulations in Fig. 6 including the temperature T , the localization lengths a and b, the charge distance d, the LO phonon energies E LO1 and E LO2 , and the local mode energy E LOC and coupling strength g/ω LOC . The local mode width ∆ and the dephasing time T 2 depend on the temperature. To reproduce the entire temperature series in Fig. 6(e) we fix the dimensions of the wave function and the local mode coupling strength for T = 300 K as given in Tab. 1. For each temperature we have to additionally determine the local mode width ∆ and the dephasing time T 2 . We find that both the ZPL and the local mode side band, which remains as narrow peak slightly above E = 2.1 eV for T = 20 K, become narrower with lower temperatures. This can already be seen when comparing the room-temperature spectra in Figs. 6(a,b) to the low-T ones in Figs. 6(c,d). In the Supplementary Material we discuss another temperature series of PL spectra taken on a different emitter, where we additionally measured the PL lifetime. In agreement with Ref. [22] we find that the life time does not depend on the temperature and that it is in the range of a few nanoseconds [17,18,19,21,27,36]. This indicates that the ZPL line width at room temperature is dominated by the pure dephasing mechanism and spectral wandering, which we combine in the fit parameter of the dephasing time T 2 . For decreasing temperatures, the ZPL width shrinks significantly as also found in Refs. [20,22]. However, the values for low temperatures between 5 K and 50 K are underestimated because the measured line width is on the order of the limited resolution of the spectrometer (100 µeV). Hence, we are not able to draw reliable conclusions about the low-temperature T 2 times from our data. The evolution of the two low-energy phonon side bands with decreasing temperature nicely supports our assumptions of their nature. The local mode results in an isolated peak at the same energy. This is not expected for phonons with continuous spectrum, i.e., for LA modes. The LA phonons lead to a low-energy side band directly attached to the ZPL, as it is well known from semiconductor QDs [55]. In the spectra in Fig. 6(e) we find that the spectral shape of the LO 1 phonon side band is additionally broadened to the high energy side compared to the simulation. When moving from high to low T , all spectral features become narrower, which makes it obvious that this broadening evolves into another maximum just below E = 2 eV. The fits for the low-temperature spectra in Figs. 6(c,d) show that our model reproduces the LO phonon side bands, despite our rough approximations for the phonon dispersion when assuming constant LO phonon energies. However, it was shown that strain strongly influences the emitters' spectra [37,28,30] and that it is highly likely that the emitters appear near the surface of the sample [23,19]. Local strain or the proximity to the surface will alter the phonon band structure and could thereby lead to variations of the spectral shape of the LO side bands. This could be a reason for the atypical broadening of the LO 1 side band for this emitter.

Photoluminescence excitation spectra
We now focus again on the distribution of ZPL energies in Fig. 2(a). We find that many of the green dots cluster around E ZPL = 2.15 eV, while the red dots cluster around 1.8 eV. These energies are approximately 150-200 meV below the respective photon energy of the laser used for excitation (dashed lines), which is the range of LO phonon energies. This suggests that LO-phonon assisted absorption could efficiently drive the emitters. Therefore, we plot a histogram of detunings between the exciting laser energy E L and the ZPL energy E ZPL in Fig. 7(a), where the colors correspond to the colors in Fig. 2(a). In total, nearly 50% of all 165 investigated emitters appear between 150 meV and 200 meV below the laser energy. The remaining emitters are almost equally distributed over the remaining 800 meV. This agrees with the findings in Ref. [26].  To confirm that the emitters are efficiently excited by LO-phonon assisted absorption, we have performed PLE measurements at room temperature to obtain information about the absorption properties of the emitters at different photon energies. The results are shown in Fig. 7(b). We chose an emitter which has a strong LO phonon sideband in the PL spectrum, as shown with the orange line. At the high-energy side of the ZPL at 2.03 eV the spectrum is again cut off by a longpass filter. We measure PL spectra for a wide range of different excitation energies from 2.09 eV to 2.3 eV as shown in the Supplementary Material. After fitting and subtracting the background and Raman lines as described in the Supplementary Material, we extract the PLE data by integrating the PL spectra either over the ZPL, i.e., from 1.95 eV to 2.05 eV (blue shaded area), or over the LO phonon sideband, i.e., from 1.75 eV to 1.9 eV (yellow shaded area). The PLE data for the ZPL and for the LO phonon sideband are shown in blue and yellow, respectively. We find that both PLE data sets exhibit a pronounced maximum between 2.2 eV and 2.3 eV. This energy range nicely agrees with the phonon sidebands of the mirrored PL spectrum (green line). This observation demonstrates that the excitation of the emitter is very efficient via the LO-phonon assisted process.
Our study shows that the LO-phonon assisted excitation of the excited defect state provides an optimized way to drive the system. Therefore we can conclude that it is conceivable that we select a specific ZPL energy ω 0 , when pumping the system with energies between ω 0 + 160 meV and ω 0 + 200 meV, i.e, via a LO-phonon assisted absorption. This provides a strategy of isolating single-photon emitters with a desired emission energy from the wide range of possible transitions.

Conclusions
We have shown that the coupling to different phonon modes plays a crucial role for localized light emitters in hBN. On the one hand, coupling to bulk phonons leads to an asymmetric broadening of the ZPL in the case of LA modes and to the appearance of prominent sidebands well separated from the ZPL in the case of LO modes. On the other hand the coupling to a local mode oscillation contributes to the broadening of the ZPL at room temperature, while the respective side band is discernible at cryogenic temperatures. By fitting measured PL spectra with our theoretical model we were able to extract parameters for the wave function geometries of the emitters. We have shown that the distance between the positive and the negative charge center may be connected to the energy of the ZPL, supporting the assignment of the ZPL energy spread with the Stark effect caused by nearby charges. It was also possible to find reasonable values for the deformation potential of electrons and holes in hBN. By measuring PLE spectra, we have demonstrated that these sidebands are also present in absorption and lead to an efficient absorption via LO-phonon assisted transitions. Finally we have shown that it is possible to preferentially select emitters with a desired transition from the wide range of emission energies possible in hBN by aiming at the LO-phonon assisted absorption. Our results lead to a deeper understanding of the fundamental properties of color centers in hBN and pave the way to a tailored control of the excitation process.

Sample preparation
Single-photon emitters (SPE) in hexagonal boron nitride are observed in hBN nanopowder (Sigma-Aldrich, grain size <150 nm). The powder is micromechanically exfoliated using scotch tape and stamped onto a Silicon substrate with a 80 nm or 270 nm thick thermal oxide layer on top. Figure 1 shows exemplary images of the nano powder sample. Panels

Photoluminescence and photoluminescence excitation spectroscopy
The photoluminescence experiments are carried out in a homebuilt confocal microscope. For excitation at 2.76 eV (450 nm) a diode laser and at 2.33 eV (532 nm) a frequencydoubled Nd:YAG laser is used. For recording the PLE spectra, the emitters are excited with a tunable laser light source based on a continuous-wave (cw) optical parametric oscillator (OPO). The turn-key system (from Hübner Photonics) covers the wavelength range of 450 nm -650 nm (2.76 eV -1.91 eV).
For the room-temperature measurements, the laser is focussed to the diffraction limit by an objective lens (numerical aperture NA=0.9), resulting in a focus size of ≈ 400 nm. The excitation power is 100 µW and kept constant for all excitation energies.
For the temperature-dependent measurements, the sample resides in a flow-cryostat cooled by liquid Helium. In this experiment, the laser is focussed on the sample using an objective lens with a NA of 0.75. In both cases, the same objective lens is used to collect the photoluminescence of the single emitters and the PL is analyzed in a spectrometer with an attached liquid nitrogen-cooled CCD camera. A longpass filter at 2.06 eV is used to remove stray light of the laser from the PL signal.

Phenomenological model
As explained in the main text we use a phenomenological fit function to extract relative spectral contributions from the phonon side bands of 165 different emitters.
The full fit function consists of three parts with where the zero phonon line (ZPL) is given by The longitudinal optical (LO) phonon side bands are assumed to have the same width as the ZPL, which leads to the function where E LO 1,2 are the LO phonon energies. Because the considered low energy (LE) phonon contributions have energies of E LE 1 = 14 meV and E LE 2 = 30 meV they are in the range of the thermal energy at room temperature k B T ≈ 25 meV. Therefore we take phonon emission (I LE j (E)) and absorption (I abs LE j (E)) processes into account. Additionally we consider two phonon processes (I 2ph LE j (E)) for these modes and combination of LO and LE processes. Therefore, the LE phonon side bands read The single LE phonon lines read where n j is the thermal occupation of LE mode j from the Bose-Einstein distribution. to reduce the number of fit parameters we assume that the weights of the two phonon processes scale like the single LE phonon processes. This leads to the fit function

Model parameters
Motivated by Ref. [1], for the LO phonon energies we consider two constant energies and for the LA phonons linear dispersions. We distinguish between the in-plane and the out-of-plane direction with The sound velocities are chosen to c z = 3.44 nm/ps [2] for the out-of-plane direction and for the in-plane direction we consider the mean value of the two high symmetry directions with c r = 16 nm/ps, which we extract from the dispersion relations in Ref. [3]. Also the dielectric constants are different in the two distinct lattice directions with [4] ε ∞ (q r ) = 4.95 , ε s (q r ) = 7.04 , ε ∞ (q z ) = 4.10 , ε s (q z ) = 5.09 We interpolate between the two given directions via ω j (q) = cos 2 (ϕ q )ω j (q r ) + sin 2 (ϕ q )ω j (q z ) ε(q) = cos 2 (ϕ q )ε(q r ) + sin 2 (ϕ q )ε(q z ) where ϕ q = atan(q z /q r ) is the angle of q with respect to the hBN-plane.

Temperature-dependent photoluminescence measurement
To confirm the conclusions from the temperature-dependent measurements drawn in the main text, in Fig. 2 we show a second data set for a different localized emitter.  Figure 2. Temperature-dependent PL spectra of one localized emitter. The temperature decreases from top to bottom (red to blue). The simulation is shown in green. The temperatures T are given next to each spectrum.
In Fig. 2(a) the measured PL spectra (color) are shown for decreasing temperatures from top to bottom. The corresponding simulations are shown as black lines. The parameters used are given in the plot next to each curve. Here only the dephasing time T 2 is adjusted for each temperature, while the other parameters are kept constant. We did not find strong indications for a contribution of the local mode, leading to a coupling strength of g = 0 in the simulations. The deformation potential couplings are the same as in the main text (D e = 0.4 eV and D h = −0.2 eV). We again find that the ZPL strongly narrows when reducing the temperature.
In addition we measured the lifetime of the PL signal for different temperatures. These measurements are carried out in the same setup as the low-temperature PL measurements (see Sec. 1.2). However, the emitters are excited by a tunable femtosecond fiber laser system (≈ 250 fs pulse length, 40 MHz repetition rate) at an energy of 2.16 eV [5]. The photoluminescence is detected using a Picoquant PDM Series single-photon counting module (timing accuracy 50 ps) and the PL decay measurement is performed with a Becker & Hickl SPC-130-EM time-correlated single-photon counting card. The excitation power is kept constant for all measurements at 150 µW. The results are shown in Fig. 2(b). We extract the lifetime T 1 by fitting the data with a convolution of the instrument response function and a single exponential decay (black curves). We find that the lifetime is almost constant T 1 ≈ 2.3 ns for all considered temperatures. These values are on the same time scale as those of defect centers in diamond [6].
In Refs. [7,8] it was shown, that the homogenous linewidth of the emitters at low temperatures is in the sub-µeV range and the corresponding T 2 in the ns range. In this case, the dephasing is dominated by the lifetime of the excited state, i.e., additional dephasing effects such as spectral wandering or pure dephasing vanish. This supports our observed trend of a strong narrowing of the ZPL with decreasing temperature. However, the values for low temperatures (5 K and 50 K) are underestimated because the measured line width is on the order of the limited resolution of the spectrometer (100 µeV). Hence, we are not able to draw reliable conclusions about the low-temperature T 2 times from our data. At room temperature, we find that the radiative lifetime is four orders of magnitude longer than the dephasing time. Therefore, we can conclude that the dephasing due to the radiative lifetime is negligible compared to other dephasing mechanisms. Figure 3 shows PL spectra of the nitrogen vacancy (NV − ) and the H3 color center in diamond in orange, respectively. The single-crystal diamond plate used for this experiment has a lateral size of 3 × 3 mm 2 and a thickness of 1 mm and is produced by a high-pressure, high-temperature process (element6) and contains approximately 200 ppm nitrogen. For the excitation of the NV − center (H3-center) we used a 532 nm laser diode (405 nm laser diode) with a power of 100 µW at the sample position. The focus diameter is 1.2 µm for 532 nm excitation and 1.4 µm for 405 nm excitation. Both spectra were measured at cryogenic temperatures. The blue curves represent the calculated spectra, where the coupling to a single local mode was considered. The model parameters with the best agreement with the experiment are listed in Tab. 2.   Fig. 3. The table shows the defect type, the ZPL energy ω 0 , the emitter dephasing time T 2 and the local mode energy ω LOC , its strength g and its width ∆.

PLE measurements
In Fig. 4(a) PL spectra for different excitation energies are shown. The ZPL maxima and the phonon sidebands reside on a background signal. For excitation energies between 2.09 eV and 2.22 eV, Raman lines appear, which shift with the laser energy. The Raman shifts of the three lines are 63 meV, 118 meV, and 168 meV. We attributed the first two to the Si substrate [9] and the third one to the LO 1 energy [10]. To extract the intensity of the ZPL and the phonon sidebands, we fit the background by an exponential multiplied by an error function to reproduce the edge of the filter. The Raman lines are modeled by Gaussians. The background-and Raman-line-corrected spectra are shown in Fig. 4(b). For small excitation energies around 2.1 eV, which are very close to the ZPL energy, The PL spectra are very weak. Therefore both PLE signals in Fig. 7 in the main text drop significantly. This surprising observation indicates that in the energetic range where excitation processes -including the local mode -should occur, the emitter is not excited. Looking at the PL spectra in Fig. 4(a) we see that for small excitation energies the ZPL overlaps with a strong Raman line, which might spoil the evaluation of the PLE data. Apart from that, the deviation between the PLE spectrum and the mirrored emission spectrum remains unclear and deserves additional investigation which is however beyond the scope of this paper.