First principles theory of the nitrogen interstitial in hBN: a plausible model for the blue emitter

Color centers in hexagonal boron nitride (hBN) have attracted considerable attention due to their remarkable optical properties enabling robust room temperature photonics and quantum optics applications in the visible spectral range. On the other hand, identification of the microscopic origin of color centers in hBN has turned out to be a great challenge that hinders in-depth theoretical characterization, on-demand fabrication, and development of integrated photonic devices. This is also true for the blue emitter, which is an irradiation damage in hBN emitting at 436 nm wavelengths with desirable properties. Here, we propose the negatively charged nitrogen split interstitial defect in hBN as a plausible microscopic model for the blue emitter. To this end, we carry out a comprehensive first principles theoretical study of the nitrogen interstitial. We carefully analyze the accuracy of first principles methods and show that the commonly used HSE hybrid exchange-correlation functional fails to describe the electronic structure of this defect. Using the generalized Koopman's theorem, we fine tune the functional and obtain a zero-phonon photoluminescence (ZPL) energy in the blue spectral range. We show that the defect exhibits high emission rate in the ZPL line and features a characteristic phonon side band that resembles the blue emitter's spectrum. Furthermore, we study the electric field dependence of the ZPL and numerically show that the defect exhibits a quadratic Stark shift for perpendicular to plane electric fields, making the emitter insensitive to electric field fluctuations in first order. Our work emphasize the need for assessing the accuracy of common first principles methods in hBN and exemplifies a workaround methodology. Furthermore, our work is a step towards understanding the structure of the blue emitter and utilizing it in photonics applications.


I. INTRODUCTION
The development of quantum communication 1-3 and quantum internet 4,5 demand the preparation and processing of tailored photonic states to carry pieces of quantum information to long distances. 6Complex photonic technologies require on-demand emission of single photons from controllable single photon emitters (SPEs).For real-life applications SPEs need to produce high purity single photons with a high emission rate, be robust against environmental disturbances, and should offer solutions for chip-scale integration in photonic devices.
Numerous material platforms, ranging from silicon through diamond to polymers, have been proposed and tested for integrated quantum photonic applications. 6Each platform offers different advantages, however, so far no material has been found that can satisfy all the needs for large scale photonic applications.Therefore, currently the most viable approach is the optimization of different processing elements in various hosts materials offering crossplatform integration capabilities. 6 the forefront of single photon emitter developments, color centers in hexagonal boron nitride (hBN) have recently gained considerable attention. 7Numerous color centers have been reported that feature sharp ZPL line at room temperature, on demand single photon emission rate in the MHz range, high spectral stability, and high quantum efficiency. 8,9ese color centers can be engineered to realize SPEs emitting at various wavelengths in a wide spectral range. 8,9Due to its layered structure, hBN also offers versatile integration capabilities for photonic applications. 7,10? -15 The latter features are attributed to a quadratic electric field dependence of the Stark shift of the blue emitter, which make the defect insensitive to weak electric field fluctuations in first order. 1513]16 In addition, the appearance of the blue emitter seems to show a close relation to another emitter in the ultraviolet (UV) spectral range 16 .8][19] This is also true for the blue emitter, as the identification of the blue emitter's microscopic structure remains elusive hindering theoretical studies and slowing down the development of related photonic applications. 13om the analysis of the blue emitter's photophysics, several important properties of the underlying atomic and electronic structure have been deduced. 15,20From the quadratic Stark effect one could infer a high, D 3h symmetry of the ground and excited state configurations. 15 addition, the stability of the center under continuous optical pumping suggests a closed shell ground state with fully occupied and completely empty defect state(s) in the lower half and the upper half of the band gap, respectively. 20Taking into considerations the conditions of the fabrication, vacancies, interstitials, antisites, and complexes derived from the combination of these defects are the most plausible candidates of the microscopic structures that can be created by irradiation in hBN. 13Related structures, such as the nitrogen interstitial and the carbon split interstitial dimer, have already been suggested in Ref. [15] as potential candidates that fulfill the above mentioned criteria.2][23][24][25] On the other hand, no comprehensive theoretical studies have been carried out for these candidates.
Here, we theoretically study the negatively charged nitrogen interstitial defect and propose it as a microscopic model for blue emitter in hBN.So far, no consensus has been achieved in the literature concerning the most favourable configuration of this defect.We show that the split interstitial configuration is the lowest energy form of the nitrogen interstitial.We carefully analyze the accuracy of popular density functional theory (DFT) methods by testing the generalized Koopman's theorem [26][27][28] for the defect states and carrying out fractional occupation number weighted electron density (FOD) 29,30 analysis on the system.We show that the excited state of the defect, exhibiting stretched nitrogen-boron bonds, cannot be described by the hybrid functional most often used for hBN.The fraction of the exact exchange contribution needs to be reduced to obtain a higher accuracy ZPL value.Furthermore, we calculate the photoluminescence (PL) spectra of the negatively charged nitrogen split interstitial and show that it resembles the experimental spectra of the blue emitter.We also study the electric field dependence of the ZPL energy and numerically demonstrate a quadratic Stark shift.We conclude that the nitrogen split interstitial defect is a promising model for the blue emitter which may facilitate advances in the fabrication and utilization of the emitter.

II. RESULTS -FIRST PRINCIPLES THEORY OF THE NITROGEN INTERSTI-TIAL IN HBN
In this section, we carry out a focused theoretical study on the nitrogen split interstitial defect in hBN.Our goal is to achieve the highest possible accuracy for the N i defect by utilizing a diverse set of state-of-the-art numerical methods.Relevant details of the methodology are incorporated into the discussions of the results in this section.The remaining technical details of the calculations are provided in the Methods section.

A. The stable atomic configuration
First, we study the stable and metastable atomic configurations of an inter-layer nitrogen interstitial in hBN.We consider the negative charge state of the system, which is expected to be the most stable for Fermi energy values close to the middle of the band gap. 22,24Various related structures have been studied in the literature before; 21-24 however, no consensus has been achieved regarding the ground state configuration.Although multiple studies identify the [0001] split interstitial as the most favored configuration, [21][22][23] a recent study found a tilted split interstitial configuration for the ground state. 24Furthermore, this tilted split interstitial configuration was proposed as a candidate for the 3.1 eV emitter. 24 order to unambiguously identify the ground state configuration of the nitrogen interstitial in hBN, we investigate a series of potentially relevant defect configurations in a multi-layered bulk systems, see Fig. 1.We first place the interstitial atom into a small orthorhombic 120-atom supercell of two layer hBN with AA ′ stacking similar to the model used in Ref. [24].Since relaxation of the defect and over-mixing of defect and host states can lead to different ground state structures 31 , we consider both semi-local PBE functional 32 and hybrid HSE functional with Fock exchange contribution of α = 0.30 and screening parameter µ set to 0.4 Å−1 (also used in Ref. [24]) to identify the ground state.We find the [0001] split interstitial, Fig. 1(a), to be favorable by 0.59 eV (0.68 eV) than the lowest energy tilted split-interstitial configuration calculated using HSE(0.3)(PBE) functional.Therefore, we could not confirm the results of Ref. [24], but our results are in-line with earlier reports of Refs [21-23]   Before discussing our results for larger, convergent supercells, we point out shortcomings of the two-layer model.To study the stability of the model, we carry out molecular dynamics (MD) simulation at 500 K using PBE functional.As a result of the atomic thermal motion, we observe sliding of adjacent BN layers from stacking AA ′ to AB1 ′ , where the B atoms are aligned on top of each other, while N atoms slide to the hollow site. 33The fully relaxed [0001]   split interstitial in the AB1 ′ stacked hBN is found to be 0.68 eV (0.95 eV) more favorable than in the AA ′ stacked hBN calculated by using HSE(0.3)(PBE) functional.This is due to reduction of the the repulsive interaction between the N atoms in the AB1 ′ stacking.From the analysis of the atomic configurations proposed in Ref. [24] and from our experiences with two-layer models, we conclude that the recently reported tilted ground state configuration 24 is possibly a consequence of unphysical layer slippage.We note here that the sliding of the BN layers results in a metastable configuration in a four-layer model.To avoid slippage of the layers in the simulations, we include at least 4 layers in our bulk models hereinafter.
According to our findings, the hexagonal 512-atom supercell comprising of four BN layers is sufficient to accurately describe the ground state configurations of the nitrogen interstitial in hBN.Using HSE(0.3)functional and the 512-atom model, the [0001] split interstitial configuration is found to be the most favorable one, see Fig. 1(a), followed by a tilted interstitial configuration with tilt angle of 17°from the [0001] axis along the hollow region and 0.51 eV higher in energy (16°and 0.50 eV with PBE), see Fig. 1(b).Another metastable interstitial configuration is obtained with 0.50 eV higher energy and pronounced tilt of 26°i n the BN bond direction (25°and 0.59 eV with PBE).Unlike the previous configurations, the 26°tilted interstitial consists of three-coordinated N atom, see Fig. 1(c).Additional calculations indicate that the energy difference between different configurations is more sensitive to the exchange fraction in comparison to separation between hBN layers and choice of van der Waals screening.
Using HSE(0.3)functional, the distance between the split interstitial N atoms and the neighboring B atoms are d 1 (interstitial nitrogen-nitrogen pair) 1. 56

B. Electronic structure
After identifying the most favourable configuration of the negatively charged nitrogen split interstitial defect in hBN, we continue our study with the analysis of the electronic structure of the defect.The two out-of-plane nitrogen atoms give rise to several defect states, some of which fall into the valence band.The ones that appear in the band gap are the fully occupied e ′′ state and the empty a ′′ 2 state, see Fig. 2. The ground state is a singlet, and the defect possesses no spin.
Optical transition between the occupied and the unoccupied defect states is allowed by perpendicular to c polarized photons.In the excited state, the nitrogen atoms form long bonds with the closest boron atoms on the symmetry axis due to the occupation of the bonding a ′′ 2 state.Consequently, the N-N distance d 1 extends by about 26% (0.42 Å ), while the B-N distance d 3 on the symmetry axis contract by 33% (0.85 Å ).We note that the excited state of the N i (-) defect is slightly distorted due to the Jahn-Teller effect that reduces the symmetry to C 2v .The three B-N distances in the middle plane (d 2 ) change from 1.58 Å, 1.58 Å, and 1.58 Å to 1.67 Å, 1.67 Å, and 1.58 Å.
Accurate first principles calculation of the optical properties of color centers in hBN has turned out to be a challenge.In order to draw reliable conclusions for the N i (-) defect, we test different methods, such as PBE and HSE(α) functionals in periodic supercell models and PBE-TDDFT and CASSCF-NEVPT2 methods on cluster models of the N i (-) defect in hBN.
In periodic 768-atom models, we use the constrained occupation method and Hellmann-Feynman forces to relax the excited state structure corresponds to e ′′ → a ′′ 2 transition.Excited states geometry is first relaxed with PBE functional and then continued with HSE(α) functional.We note that relaxation of the excited state with HSE(α) (α ≈ 0.25 -0.35) is cumbersome due to the change of ordering of the defect states and related structural instabilities.In our 261-atom cluster model of 5 adjacent hBN flakes, see Fig. 3, we initially optimize the ground state geometry at PBE level of theory, using def2-SVP basis set.During the optimization, owing to its computational demands, only the middle 13 atoms are relaxed in each layer plus the additional interstitial nitrogen.On the obtained equilibrium geometry, we calculate the vertical excitation spectrum with PBE-TDDFT and CASSCF-NEVPT2 methods, requesting 3 roots (ground state and 2 degenerate excited states in all cases).Then, we relax the structure of the excited state by PBE-TDDFT, following the energy gradient of the corresponding root.The electronic energy of the excited state is calculated from a second vertical excitation spectrum calculation at the relaxed geometry.Using a larger basis set would have been too expensive, however to obtain more accurate ZPL, def2-TZVPD basis is applied on the middle 4 atoms of each layer plus the interstitial nitrogen while def2-SVP basis set on the remaining atoms.More details on the excited state calculations and the models can be found in the Methods section.
The adiabatic energy differences (E AD = E ES − E GS ) calculated for the various models are provided in Table I.The true ZPL energy includes contribution from zero point energies of the phonon modes in the ground and excited states ( As can be seen, the adiabatic energy differences ranges in a wide energy interval starting at 2.38 eV up to 3.20 eV.The wavefunction-based NEVPT2 method predicts the highest energy difference of 3.20 eV.As demonstrated for other systems [34][35][36] , the CASSCF-NEVPT2 method can provide accurate results when excited states of multi-reference character are considered.However, for Translucent green denotes the FOD surface at 0.005 e/Bohr 3 and T=5000 K. the N i (-) defect, both the ground state and the lowest energy singlet excited state possess a clear single-reference character.The leading electron configurations visualized in Fig. 2 have 100.0%and 98.9% contribution in the ground and excited state wavefunctions, respectively, if 4 electrons and the 3 frontier orbitals form the active space.We note that similar wavefunction characters are obtained for systematically extended active spaces.Therefore, the accuracy of the NEVPT2 method in this case is comparable to the MP2 method 37 , which is expected to be outperformed by a carefully chosen single-determinant DFT method.
As an alternative investigation of the reliability of single-reference methods, we adopt the fractional occupation number weighted electron density (FOD) analysis method 29 from quantum chemistry.Using this method one can identify potential static correlation features of the obtained DFT ground states, which might call for more sophisticated multi-reference methods.To this end, we carry out calculation on the same 261 atom cluster model, Fig. 3(a), using PBE functional and def2-SVP basis set.
The FOD analysis performed on the ground state geometry reveals no FOD density, see Fig. 3(b).For the excited state geometry, however, some FOD density is predicted, which is found to be localized around the defect center on the s N-B bonds stretched along the symmetry axis, see Fig. 3(c).The shape of the FOD surface resembles the union of the two HOMOs and the LUMO with no detectable contribution from other orbitals.Thus, static correlation is only expected to be present in the aforementioned 4-electron 3-orbital (4,3) active space.As CASSCF(4,3) calculation did not result in mixed states, it can be concluded that the obtained FOD plot indicates merely the quasi-degeneracy of the frontier orbitals and does not stem from actual multireference character.Accordingly, we continue our analysis using TDDFT and constrained-DFT methods.
The PBE-TDDFT result obtained in a 261-atom cluster using def2-SVP basis set and the PBE constrained DFT result obtained in a 768-atom periodic model are provided in Table I.
As can be seen, the two methods agree very well with each other indicating a proper basis set selection for the cluster model.On the other hand, PBE-DFT and the NEVPT2 results deviates by ≈0.8 eV.Since semi-local exchange-correlation functionals tend to underestimate and the NEVPT2 (∼MP2) tends to overestimate excitation energies, we set our focus to hybrid DFT methods, which can be grasped as an interpolation between (semi-)local density functional theory and wavefunction theory.
In order to further narrow down the set of applicable methods, we select the HSE06 38 hybrid-DFT functional for further investigations; nevertheless, we varied the mixing parameter (α, indicating the proportion of HF exchange in the functional formula) and thereby generated different HSE(α) methods.As the next step, we calculated the ZPL energy at HSE(α) levels of theory in 768-atom periodic supercell model.
The most commonly used mixing parameter α range between 0.3 and 0.35 in hBN.These values are set by adjusting the theoretical band gap to the experimental one. 22While this method improves the description of host states and localized orbitals of intrinsic defects, it may not work for all defects, especially when stretched bonds are observed, see Fig. 2(b).
In order to test the accuracy of the functional for the description of a given state, we used the generalized Koopman's theorem 26-28 (ionization potential theory 39-41 in other terminology).According to this theorem, the eigenenergy of the highest occupied Kohn-Sham orbital and the ionization energy of the system should be equal when the exact exchange-correlation functional is used.Since the exact DFT functional is not known, approximate functionals may fail to fulfill this theorem.This failure manifests itself as a difference between the KS eigenenergy and the ionization energy.The non-Koopman's energy E NK measures this as where E N +1 and E N is the total energy of the system consisting of N +1 and N electrons, and ε HO is the Kohn-Sham energy of the highest occupied orbital in the N + 1-electron system.
In practice, the non-Koopman's energy calculated in periodic boundary conditions depends also on technical parameters of the utilized model that needs to be taken into account, where E ′ NK is the non-Koopman's energy calculated in a periodic supercell model and ϵ corr is a correction term that cancels finite-size effect of the energy term originating mainly from spurious electrostatic interaction in periodic boundary conditions.Since the initial N + 1 electron system is the negative charge state of the N i defect, ϵ corr incorporates the charge correction of both ε HO and E N +1 .Both E ′ NK and the ϵ corr energy terms depend on the perpendicular and parallel to c axis extension of the supercell, L ⊥ and L ∥ respectively, as well as the charge state q of the defect.
The non-Koopman's energy defined in Eqs.(1) and (2) is a good indicator of the accuracy of DFT functionals for localized defect states.When a functional with tuneable parameters is used, such as the HSE(α) functional, adjustment of the free parameter may be utilized to decrease the non-Koopman's energy and thus to improve the description of the defect orbitals.This strategy has been successfully employed for several other point defect systems, see for instance Ref. [28].On the other hand, we must first calculate the ϵ corr term, which may be comparable in absolute value with the non-Koopman's energy itself.To proceed with the analysis on the adiabatic excited state-ground state energy difference of the N i (-) defect, we carefully investigate the value of the ϵ corr charge correction term.

Charge correction of the N i (−) defect in hBN
Charge correction of layers and surfaces is an active field of research today, however, the literature on charge correction methods of bulk models is fairly established. 42,43The most accurate, but computationally expensive approach is to carry out a finite-size scaling test and extrapolate energy terms to the single defect limit.To quantify ϵ corr in Eq. (2), we apply this procedure and consider 11 supercells of different sizes, each of which includes a single N i defect either in the negative (N + 1 electron system) or in the neutral charge state (N electron system).Due to the large computational demand of finite-size scaling tests, we carry out this study using PBE functional.We further assume that the results of finite-size scaling are independent of the choice of the exchange-correlation functional.Interestingly, the scaling of the Kohn-Sham energy and the total energy shows an unexpected behavior as depicted in Fig. 4. Considering the L ⊥ perpendicular to c lateral size dependence of both the Kohn-Sham energy and the ionization energy, we observe converging curves that take finite values at the L ⊥ → ∞ limit.This behaviour is typical for bulk models.On the other hand, for the L ∥ parallel to c lateral size of the supercell, we obtain a linear dependence, see Fig. 4, that diverges at the L ∥ → ∞ limit.This odd behaviour is unexpected for bulk model of several layers of hBN sheets in periodic boundary conditions.Similar linear dependence is, however, observed in single layer hBN sheet in vacuum in periodic boundary conditions 44,45 , where the electrostatic interaction of the infinite charged layers diverges.Therefore, we conclude that the charge correction of Kohn-Sham and total energy terms in bulk hBN model cannot be described by charge correction methods developed for conventional bulk semiconductors, such as diamond and silicon.The charge correction should account for the electrostatics of a localized charge in a 2D layer immersed in the media of relative permittivity other than 1 (due to the presence of other pristine hBN layers).
In order to model the systems, we express finite-size and charge dependent energy terms where E refers to either a Kohn-sham energy or the total energy of the system, E 0 is the finite-size effect-free energy value that corresponds to the limit of a single point defect in the material, a, b, and c are the free parameters that we will use to fit Eq. (3) to the calculated energy curves.S is the surface of the hBN layers in our finite periodic models.Note that for infinite layers, i.e. S → ∞, the last term on the right hand side of Eq. (3) vanishes.
For finite supercells the last term is of high relevance, though.Since we are interested in the non-Koopman's energy, we calculate the supercell size dependence of this energy term in the ground state and fit Eq. (3) directly to the obtained curve, see Fig. 5.
Considering the perpendicular to c lateral size dependence of the non-Koopman's energy at L ∥ = 26.2Å, we observe converging tendency, see Fig. 5(a), which is best fitted by a 1/L 3 ⊥ function using points beyond 15 Å.The smallest supercell is apparently an outlier of the general trend expectedly due to the overlap of the defect states.As can be seen in Fig. 5(a), the non-Koopman's energy only slightly changes beyond L ⊥ = 20 Å, thus we consider supercells of this later size to be convergent in the perpendicular to c direction.Therefore, for L ⊥ ≥ 20 Å we set parameters a and b to zero in Eq. (3), which approximation causes an error no larger than 15 meV.Parallel to c supercell size dependence of the non-Koopman's energy for L ⊥ = 20.0Å is depicted in Fig. 5(b).The points can be perfectly fitted by a linear curve, which intersects the energy axis at 0.145 eV for L ∥ = 0. From the extrapolation to zero, we can obtain the finite-size effect free non-Koompan's energy for the considered charged layered system, see Ref. energy, see Eq. (2), is ε corr = E NK − E ′ NK (−1, 20.0, 20.0) = 0.134 eV.The periodic model related electrostatic interaction between the negatively charged defect and the homogeneous positive background depends negligibly on the small variation of the defect states due to the choice of the functional.Therefore, we use the same correction for different HSE(α) functionals in the same supercell of 768 atoms.

Optimization of the HSE(α) functional for the electronic structure of the N i (-) defect
Having the charge correction of the non-Koopman's energy determined, we compute the finite-size effect free non-Koopman's energy value according to Eq. (2).In order to evaluate the accuracy of HSE(α) functional, we calculate the mixing parameter α dependence of the non-Koopman's energy for the ground and the excited states with and without structural relaxation, see Fig. 6.When no relaxation is taken into account, we use optimized ground and excited state structures as obtained by the HSE(0.30)functional.
First, we study the accuracy of the HSE(α) functional on the ground state of the negatively charged N i defect.As can be seen in Fig. 6, E NK = 0 is achieved at α = 0.285.
Note that geometry relaxation has only negligible effects on the non-Koopman's energy, see  To this end, we calculate the mixing parameter dependence of the adiabatic energy difference of the ground and excited states, see Fig. 7 and Table I.As can be seen, the energy difference decreases as the exact exchange contribution is reduced.At α = 0.285 and α = 0.132 the ground state and the excited state is described with high accuracy, respectively.The adiabatic energy differences corresponding to these mixing parameter values are 3.27 eV and 2.84 eV.We assume that the true value of the adiabatic energy difference could be found in this interval.On the other hand, due to the different behaviour of the ground and excited states, the accurate value of the adiabatic energy difference cannot be obtained by a HSE(α) functional using a single value of the mixing parameter.At the same time energy differences cannot be computed using two different funtionals.Therefore, we define the most plausible value of the adiabatic energy difference by the mean of the values obtained with the ground and excited state optimized functionals, which is E AD = 3.06 eV.
Due to the approximately linear dependence of the adiabatic energy difference on the mixing parameter, the same value can be obtained by HSE(α), where α is the mean value of the optimized mixing parameters.The expected error margin of the adiabatic energy difference layers, C is the size of the layers, and D refers to the region which was optimized in the layers.For example, 5-3-S-M means that a 5 layer model with the smaller layer was utilized and only 13 atoms per layer was optimized in the 3 middle layers plus the interstitial N.
Convergence of the photoluminescence spectra as a function of the cluster size and relaxation strategy is depicted in Fig. 8a.Due to the relatively large change of the position of the nitrogen atoms and the first neighbor atoms along the symmetry axis, see Fig. 2, the vertical gradient approximation 48,49 is utilized.In case of the smallest model of three layers (3-3-S-S), the intensity of the second peak is larger compared to the first peak associated to the ZPL.Similarly, the 3 layer model with the M region optimized (3-3-S-M) exhibits a large sideband.On the contrary, the two spectra from the largest 5 layer models (5-3-S-M and 5-5-S-M) show a strong ZPL emission and much less pronounced sideband.The same can be said if the larger layer is used (L) which implies that convergence of the PL spectrum can be reached in the flake models.Note that in Fig. 8a all transition lines, including the ZPL, were broadened by an empirical 300 cm −1 linewidth.
For the convergent model (5-5-S-M) we obtain a Huang-Rhys factor of 0.491, and radiative lifetime of 54 ns using a refractive index equal to 2.13.The corresponding transition dipole moment vector is (0.37574,-0.19576,-0.00215) in atomic unit, where the z direction is parallel Note that these interesting modes are localized around the defect center.
to the c axis of hBN.It must be noted that due to the large difference between the ground and excited state, our investigation is limited to using the vertical gradient approximation that only uses the ground state geometry and Hessian.The transition dipole moment between the two states could differ considerably from what we obtain from PBE-TDDFT on the excited state geometry.In principle, we expect the lifetime to be overestimated in our calculations.
From the analysis of the partial Huang-Rhys factors, we deduce that the dominant local vibrations responsible for the phonon sideband are i) the movement of the atoms surrounding the defect causing the Jahn-Teller distortion, ii) the stretching and contraction of the N-N bond in the defect, and iii) the movement of the B atoms from the neighbouring layers towards the defect, see Fig. 9.These vibrations account for the geometry difference between the ground and excited states.
Lastly, we calculate the zero-point energy (ZPE) contribution to the adiabatic energy difference obtained in the previous section.From the sum of the two contributions we compute the ZPL energy.The ZPE contribution is approximated by numerical frequency calculations where a 5 layer model with 261 atom is used and the middle 13 atoms are relaxed in each layer plus the interstitial nitrogen.For the ZPE contribution we obtain -0.10 eV meaning that the zero point energy contribution of the local vibrational modes is larger in the ground state than in the excited state.This is understandable as stretched bonds are formed in the excited state configuration that give rise a softer potential and lower vibrational frequencies.We assume in this paper that the zero point energy contribution does not significantly depend on the electronic structure calculation method used, thus we use the same zero point energy contribution in all ZPL value approximations.The third row of Table I provides the ZPE corrected adiabatic energy differences as obtained with different methods.

D. Electric field dependence of the ZPL energy
Finally, we numerically investigate the electric field dependence of the of the ZPL line of the negatively charged nitrogen split interstitial in hBN.As discussed in Ref. [15 and   20], defects exhibiting high D 3h symmetry in both the ground and excited state may be less sensitive to external electric fields and only second or higher order electric field dependence is expected.The negatively charged N i defect does exhibit D 3h symmetry in the ground state, as also discussed in Ref. [15]; however, the excited state is Jahn-Teller unstable that lowers the symmetry to C 2v .Therefore, it is a question whether the desirable quadratic dependence of the ZPL energy on the electric field is preserved even in such a low symmetry.
Our preceding study 15 suggests that the quadratic term dominates the electric field dependence of the ZPL, however, due to the applied large-scale model and consequent numerical uncertainties the results need to be reconsidered by an improved independent model.Here, we apply a cluster model, atomic basis functions, and ORCA suit opposed to the slab model, plane wave basis set, and VASP package utilized previously 15 .
To study the effect of the electric field on the zero-phonon line of the N i (-) defect in hBN, a static field is applied perpendicular and parallel to the layers of the cluster model with varying strength ranging from -0.05 to 0.05 V/ Å.For the calculations, we use PBE- TDDFT method, which somewhat underestimates the ZPL energy, see Table I, however, it can capture to the polarization of the defect orbitals induced by the external electric field.
We expect our results to be qualitatively accurate.
As can be seen in Fig. 10, perpendicular to the layers we find that the ZPL energy decreases quadratically with increasing strength of electric field.However, if the electric field is parallel to the layers, the ZPL energy depends almost linearly on the strength of the electric field.This behaviour is expected from the changes in dipole moment (∆µ = (−0.67,0.42, −0.01) Debye) and polarizibility (∆α) in both directions.In the parallel to layer direction we obtain 55 Å3 for ∆α, and 0.01 Debye for ∆µ, while in perpendicular direction 17 Å3 and 0.67 Debye, respectively.We do not find a perfectly linear curve in x direction due to the non-zero quadratic contribution to the Stark shift.

III. DISCUSSION -PLAUSIBLE MODEL FOR THE BLUE EMITTER IN HBN
Finally, we compare our numerical results on the negatively charged nitrogen [0001] split interstitial to the experimental reports of the blue emitter in hBN 20 .In summary, we propose the N i (-) defect as a working model for the blue emitter, although not all the experimental observation can be explained with the current model.Below, we list all the known pros and cons of this assignment.

Pros
The blue emitter in hBN can be intentionally created by electron irradiation.As a primary damage of high energy electron bombarding, boron and nitrogen vacancies and interstitials are created.Complex defect formation is expected as a result of secondary processes, such as migration of the created defects and possible recombination with other defects in the lattice.The nitrogen split interstitial is a primary radiation damage and it is expected to be formed in substantial concentrations.
Regarding the optical properties of the blue emitter, the underlying microscopic structure is expected to possess D 3h symmetry, while the relevant electronic structure of the defect should consists of a fully occupied defect state in the lower half of the band gap and an unoccupied defect state in the upper half of the band gap.The latter conditions are deduced from the stability of the defect under continuous excitation.The negatively charged nitrogen split interstitial defect possesses a fully occupied e ′′ state and an empty a ′′ 2 state inside the band gap giving rise to a singlet ground state.This electronic structure is in agreement with the expectations derived from the observed properties of the blue emitter.
The blue emitter is known to emit light mostly along out of plane directions with inplane electric field polarization.Interestingly, the perpendicular to c polarization of the blue emitter is not homogeneous, there is a preferential in-plane polarisation direction. 15is observation also implies that the three-fold rotation symmetry of the defect should be violated either in the ground or the excited state, which has been overlooked previously.
To explain this behaviour, we draw attention to the Jahn-Teller instability of the excited state of the N i (-) defect.Since an electron is excited from the e ′′ state to the a ′′ 2 state, the symmetry is reduced to C 2v in the excited state due to the Jahn-Teller distortion that gives rise to a clear single axis polarization pattern in the calculated transition dipole moment.This is in agreement with the observations.The blue emitter exhibits a well resolved ZPL at 2.844 eV (436 nm) 15 .While it is challenging to accurately calculate the ZPL energy of the nitrogen interstitial, our prediction agrees within the error margin of our calculations with the ZPL energy of the blue emitter.
We obtained 2.96 eV (419 nm), which is 12 meV (17 nm) blue shifted compared to the blue emitter's ZPL.This deviation is well within the estimated error bar of ±0.2 eV of our calculations.
Finally, we discuss the PL spectra of the blue emitter and the N i (-) defect.As can be seen in Fig. 8, the experimental and the convergent theoretical PL spectra agree very well with each other.Both spectra exhibit dominant emission in the ZPL and a vanishing sideband.
For the N i defect we obtain a Huang-Rhys factor of 0.491 and Debye-Waller factor of 0.61.
Due to the Jahn-Teller distorted excited state, the N i (-) defect possesses a non-zero electric dipole in this state.The dipole is in the plane of the hosting hBN layer.Depending on angle δ of in-plane electric field and the electric dipole of the defect, one may observe either linear (δ = {0 • , 180 • }), or quadratic (δ = {90 • , 270 • }), or mixed electric field dependence.This observation may explain the differences of the experimental results reported in Refs.[15,   50, and 51].

Cons
For perpendicular to plane electric fields we found quadratic Stark shift, while the paper of Zhigulin et al in Ref. [15] reports on a weak linear dependence.We attribute this difference to perpendicular to plane disturbances, for instance due to bending of the hBN sample.Such a distortion can lift the in-plane reflection symmetry of the system, which in turn enhances an out-of-plane electric dipole component.
The blue emitter possesses an excited state lifetime as short as ∼ 2.15 ns. 20Our computational results suggest a longer radiative lifetime in the range of 54 ns.We note that the calculations account only for radiative processes, thus the experimental excited state is expectedly shorter than the theoretical estimates.However, the deviation between the two values cannot be only explained by nonradiative processes.We anticipate that the deviation is due to the approximations used in the calculations.
Furthermore, according to recent study on the fabrication of the blue emitter, there seems to be a correlation between the appearance of the blue emitter and the UV emitter in hBN. 15,16,20The latter has been assigned to the the carbon dimer recently. 17Since the nitrogen interstitial and the carbon dimer structures have no common roots, the connection between the two defects is not obvious.We anticipate that not the formation of the atomic structures are related, but instead the stability of the optical bright state is affected similar ways.For instance, due to their similar electronic structure, the position of the Fermi level may affect the two defects similarly.This may explain the presumed relationship of the appearance of the two defects.

IV. METHODS
A. Further details on the calculation of periodic models using VASP The density functional theory calculations for bulk models are performed using Vienna Ab-initio simulation package (VASP) 52 , within the projector-augmented wave (PAW) method 53 with a plane-wave cutoff of 500 eV.The exchange correlation is described by the generalized gradient approximation of Perdew, Burke, and Ernzerhof (PBE) 47 , and we also consider the screened hybrid functional of Heyd-Scuseria-Ernzerhof (HSE06) 38 with modified exact exchange fraction of 0.1-0.35 and screening parameter of 0.4.The van der Waals interaction is included using Grimme-D3 method with Becke-Johnson damping 54,55 .
The unit cell of bulk hBN is optimized with 15×15×5 k-point grid.We obtain lattice parameters a = 2.50 Å, c = 6.56 Å, which are in agreement with experimental values of a = 2.50 Å, c = 6.60 Å measured at 10 K 56 .In the case of defect stability study, we consider orthorhombic supercell of 120 atoms (5a, 3a+6b, c) as well as hexagonal supercell containing 512 or 768 atoms.The reciprocal space is sampled with Γ-point.Atom relaxation is carried out until the forces are less than 0.01 eV/ Å (0.025 eV/ Å) with PBE (HSE(α)) functional.
We also include the non-spherical contributions within the PAW spheres.
B. Further details on the numerical studies of cluster models using ORCA The DFT calculations of the cluster models are made by ORCA 5.0.3 57 .Most of the calculations are made with PBE functional and def2-SVP basis set, if otherwise not noted.During all calculation with ORCA resolution of identity (RI) approximation are used to the Coulomb term and the fitting basis are applied through the AutoAux feature.Geometry optimizations are made with default settings.For excited states different number of roots are required depending on the setups, and the corresponding root is chosen to be followed which represents the correct excited state.FOD analysis is performed at 5000 K as suggested in the seminal article 29 .Numerical frequency calculations are performed on only those atoms which were originally optimized at the particular model setup with steps of 0.005 Å3 .The PL spectra are determined based on the Hessians obtained from the numerical frequency calculations, using the vertical gradient approximation 48,49 and a 300 cm − 1 linewidth.

APPENDIX A
In order to demonstrate that our cluster model is capable of providing convergent predictions for the studied system, we calculate the PBE-TDDFT ZPL energy with varying flake size, number of stacked layers, and geometrical optimization scheme.The corresponding numerical results are collected in Table II.Note that ZPL energies in this study do not include ZPE contribution.
Interestingly, we find that the ground state is relatively prone to the relaxation, the energy gap increases by granting more freedom for geometry relaxation than the minimalist scheme S. By adding external layers to the minimal model of three stacked layers, the in-layer planarity of the excited state is preserved yielding some 0.10-0.2eV energy gain.This stabilization effect is found not to be sensitive to the relaxation of the additional covering layers.In conclusion, independently of the investigated cluster setups, the lowest singlet excitation is consistently found at around 2-2.5 eV.Based on the ZPL results in the following analysis we study the 5 layer model with smaller layers and optimize all of them with the M strategy.The used def2-SVP is a minimal basis set, but using larger basis set would be too expensive.To address this issue, we also study how ZPL changes if larger basis sets are used on the center atoms.In our convergence test, def2-SVPD, def2-TZVP, and def2-TZVPD is applied only on i) the two N-s and two B-s in the center (4 atoms), on ii) 4 atoms at the center of the 3 middle layer plus the interstitial N (13 atoms) or on iii) 4 atoms of the 5 layer plus the interstitial N (21 atoms).The ZPLs are collected in Table III.Using diffuse functions causes less than 0.1 eV change in the ZPL, however the application of TZ quality basis set increases the ZPL from 2.20 eV to 2.28-2.38 eV.It can be also observed that including the outer layers have negligible impact.As a convergent PBE-TDDFT adiabatic ZPL value for the cluster model we obtain 2.38 eV in good agreement with the results of constrained-DFT calculations in a periodic model, see Table I While the model describes an additional interstitial N atom, each small (S) or large (L) layer is defined by 52 and 94 atoms, respectively.In the small (S), medium (M), large (L) relaxation scheme besides the hydrogens the central 4, 13, 37 second-row atoms per layer are taken into account as illustrated in Fig. 3a, respectively.Additionally, ZPL-s with HSE(0.We acknowledge KIF Ü for awarding us access to computational resource based in Hungary. FIG. 1. Atomic structure of negatively charged nitrogen interstitial defect in hBN.(a) The split interstitial with interlayer bonding is the energetically favored configuration.We note that there is no covalent bond between the two nitrogen atoms of the split interstitial defect.(b) and (c)denote the metastable configurations with intralayer bonding, and formation energies 0.51 eV and 0.50 eV (0.50 eV and 0.59 eV) calculated using HSE(0.3)(PBE) higher than the split interstitial, respectively.Nitrogen (boron) atoms are shown as blue (pink) spheres.

FIG. 2 .
FIG. 2. Kohn-Sham energy levels of the ground and excited states.LUMO, HOMO and HOMO-1 frontier orbital energies are highlighted in red, green and blue, respectively.The shape of the corresponding frontier orbitals are also plotted in both side and top view.Note that for better visibility only the central three layers of the molecules are shown.
FIG. 3. a, shows the structure of the larger layer from top view.Note that the second-row atoms included in the smaller layer are enclosed by the black hexagon while the grey, pink and blue balls denote the position of the hydrogen, boron and nitrogen atoms, respectively.The translucent red, yellow, and green highlighting represents the minimal and the additional subsets of atoms included in molecular geometry relaxation, respectively.For more details see main text.FOD analysis for the ground state and the excited state geometry presented in b and c subplot, respectively.
FIG. 4. Finite-size scaling of the Kohn-Sham energy and the total energy difference.∆ε KS is the Kohn-Sham energy of the highest occupied defect orbital measured from the energy of the valance band maximum of the pristine supercell and ∆E tot = E N +1 − E N is the ionization energy of the negatively charged N i defect.Dark blue line with circles and dark red line with diamonds show the scaling of the energies as a function of the in-plane supercell size L ⊥ , while light blue line with squares and pink line with triangles show the energies as a function of the perpendicular to plane supercell size L ∥ .

FIG. 5 .
FIG. 5. Determination of the charge correction of the Non-Koopmans energy.a and b depict the in-plane (L ⊥ ) and the parallel to c (L ∥ ) supercell size dependence of the Non-Koopman's energy (E ′ N K ).ϵ corr is the charge correction of the Non-Koopman's energy at L ∥ = 20 Å.

Fig. 6 .FIG. 7 .
Fig.6.The obtained optimal mixing parameter is close to α = 0.3 used in the analysis of the ground state properties, thus our results presented for the stability of the N i defect are valid.On the other hand, for the excited state we observe a different behaviour.At α > 0.3 the finite-size effect corrected non-Koopman's energy takes a significant negative value, indicating that the frequently used HSE({0.3,0.32,0.35})functionals are not appropriate for the description of the excited state of the N i (-) defect.For the excited state, the non-Koopman's energy vanishes at α = 0.132 (α ≈ 0.154) when structural relaxation is taken (not taken) into account.Here, structural relaxation makes a difference in the non-Koopman's energy in contrast to the case of the ground state.These results clearly indicate the need for the reduction of the exact exchange contribution in the excited state.Since accurate description of the ground state and the excited state requires two different functionals, HSE(0.132) and HSE(0.285), the adiabatic energy differences of the states cannot be calculated with the highest accuracy.On the other hand, we can use

FIG. 8 .
FIG. 8. Left: calculated PL spectra with different geometry optimization setups.Right: The best calculated PL spectrum compared to experimental results 20 .Partial Huang-Rhys factors of each mode are also presented with bars.Linewidth is set to 300 cm −1 .For better comparison of the theoretical and the experimental PL spectra, the theoretical spectrum is shifted to match the experimental ZPL.

FIG. 9 .
FIG. 9. Visualization of the four most dominant phonon modes based on the Huang-Rhys factors.

FIG. 10 .
FIG. 10.The electric field dependence of the ZPL.The ZPL energy is calculated on PBE-TDDFT level of theory in a 261-atom cluster model.Due to the approximations used, the absolute value of the ZPL energy is underestimated.
knowledges the Australian Research Council (CE200100010, FT220100053) and the Office of Naval Research Global (N62909-22-1-2028) for the financial support.V.I. also appreciates support from the Knut and Alice Wallenberg Foundation through WBSQD2 project (Grant No. 2018.0071).The calculations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at the National Supercomputer Centre (NSC).

TABLE I .
Adiabatic energy differences and ZPL values as obtained with different computational methods.The values in parentheses are given in nm, while all the other values are in eV.The zero point energy contribution (ZPE) contribution to the adiabatic energy difference is -0.10 eV.Bold font indicates our most accurate theoretical model, adiabatic energy difference, and ZPL.

TABLE II .
. The effect of different number and size of layers, and optimized regions on the TDDFT-PBE and single point HSE ZPL values.ZPL energies in this table do not include ZPE contribution.

TABLE III .
25)functional are determined on the geometries obtained with PBE.These results are consistently higher around 1 eV regardless The effect of larger basis sets on selected atoms on the PBE-TDDFT ZPL.