Ultraviolet Quantum Emitters in Hexagonal Boron Nitride from Carbon Clusters.

Ultraviolet (UV) quantum emitters in hexagonal boron nitride (hBN) have generated considerable interest due to their outstanding optical response. Recent experiments have identified a carbon impurity as a possible source of UV single-photon emission. Here, on the basis of first-principles calculations, we systematically evaluate the ability of substitutional carbon defects to develop the UV color centers in hBN. Of 17 defect configurations under consideration, we particularly emphasize the carbon ring defect (6C), for which the calculated zero-phonon line agrees well the experimental 4.1 eV emission signal. We also compare the optical properties of 6C with those of other relevant defects, thereby outlining the key differences in the emission mechanism. Our findings provide new insights into the strong response of this color center to external perturbations and pave the way to a robust identification of the particular carbon substitutional defects by spectroscopic methods.

Single-point defects in two-dimensional (2D) hexagonal boron nitride (hBN) play a vital role in the optical properties of the host and hold great promise for quantum information technologies and integrated quantum nanophotonics.[1−8] In particular, color centers in hBN are responsible for ultrabright single-photon emission at room temperature with a wide range of emission wavelengths.[8,9] Recent experiments demonstrated the versatile properties of the defect emitters in 2D hBN, such as strain- and electric field-dependent emission,[5,10−12] high stability under high pressure and temperature,[13−15] and initialization and readout of a spin state through optical pumping.[2,16] Other studies have shown a successful engineering and coherent control of a single spin in hBN,[3] while room-temperature initialization and readout have also been realized.[2,16] Of several photoluminescence (PL) signals from the color centers in hBN, multicolor single-photon emissions have been detected around 1.6–2.2 eV. Dozens of studies have been performed to determine the possible origin based on simple defect configurations.[1−5,8,17,18] In addition, a strong ultraviolet (UV) emission at close to ∼4.1 eV has received a great deal of attention.[19−23] The single-photon emission associated with these bands indicates that it should originate from a point defect.[6,24] However, despite various attempts, the atomistic origin of the UV emission in hBN is still under debate. In particular, due to the similarities with the carbon-doped hBN samples (mostly due to the PL lifetime of ∼1.1 ns[19,25]), carbon is thought to contribute to the formation of the PL signal.[21,26] Despite the fact that some of the proposed configurations exhibit excitation energies around 4 eV,[17,21,23,27−30] many of their key properties, including the stability, electronic configuration, and vibronic properties, were not considered. Recently, additional lines were observed in the range of 4.1–4.2 eV and isotopically controlled carbon doping is employed to determine the role of the carbon impurity in the 4.1 eV emission.[23] In particular, the additional lines, distinct from the previous 4.1 eV emission, show strong PL intensity with a clear temperature dependency.[15] These findings motivated us to carry out a systematic theoretical study to reveal the role of substitutional carbon defects in the formation of the UV single-photon emitters in hBN.

In this paper, we analyze 17 configurations of substitutional carbon defects and systematically address their thermodynamic properties. Among those, we identify a six-carbon ring defect, in which the carbon atoms substitute one BN honeycomb of hBN lattice, as one stable defect configuration. It is noteworthy that this defect has already been unambiguously identified via annular dark field scanning transmission electron microscopy (ADF-STEM)[31,32] and can be intentionally introduced into the lattice with atomic precision by the focused electron beam.[32] We show that this color center emits light due to strong electron coupling with E-phonon modes, caused by the product Jahn–Teller effect. More specifically, the respective symmetry lowering is found to activate a forbidden transition through an intensity borrowing mechanism from a higher-lying bright state. We further calculate the zero-phonon line (ZPL) energy, luminescence spectrum, and radiative lifetime and found them to be in excellent agreement with the experimental observations for the 4.1 eV emission. In addition, we discuss the possibilities of distinguishing between different carbon configurations on the basis of the 13C isotopic shift in ZPL and sideband and on the basis of a different response to the applied strain.

The calculations were performed on the basis of the spin-polarized density functional theory (DFT) within the Kohn–Sham scheme as implemented in the Vienna ab initio simulation package (VASP).[33,34] A standard projector-augmented wave (PAW) formalism[35,36] was applied to accurately describe the spin density of valence electrons close to nuclei. The screened hybrid density functional of Heyd, Scuseria, and Ernzerhof (HSE)[37] was used to optimize the structure and calculate the electronic properties.

The calculations with the second-order approximate coupled cluster singles and doubles model (CC2)[38] and the algebraic diagrammatic construction method [ADC(2)][39] were performed with the Turbomole code.[40,41] The results of time-dependent (TD) DFT and n-electron valence state perturbation theory [NEVPT2(4,4)][42] were obtained with the ORCA code.[43] We used the cc-pVDZ basis set[44] and considered the PBE0 density functional[45] for TDDFT. The periodic TD-PBE0 calculations were performed with Quantum Espresso.[46] The details of modeling, calculation parameters, and the computation of formation energies together with the charge correction[47] are discussed in Supplementary Note 1.

First, we systematically analyzed the thermodynamic properties of the carbon defects in hBN. Because the experimental PL signal features a short radiative lifetime, we focused on only those arrangements in which the carbon atoms are closely packed within a single honeycomb. The delocalization of defect orbitals should naturally decrease the excitation energy (see Figure 2 of the Supporting Information); therefore, larger defect complexes were not considered. The resulting structures of 17 distinct C configurations are shown in Figure a. For those, we evaluated the formation energy diagrams and charge transition levels (CTLs), which are plotted in Figure b and Figure 1 of the Supporting Information. The formation energies for the defects with an unequal amount of substituted B and N can be largely decreased by selecting the appropriate growth conditions. However, for a given number of carbon atoms, we always observed that the most stable configurations represent the confined C clusters, where the carbon atoms are arranged in a continuous chain. Importantly, to prevent a photoionization process, a UV quantum emitter should maintain a stable charge state. This condition is observed for the defects with an even number of carbon atoms (namely, CNCB, 2CNB, 4Cchain, 4Cpair, and 6C); they possess a highly stable neutral charge state over the energy range, exceeding the ionization threshold. By contrast, the defects with an odd number of carbons rapidly change their charge states across the formation energy diagrams because of their radical nature. Our calculations provide a low formation energy of 2.17 eV for carbon dimer CNCB, which is quite consistent with the previous reports.[27,48] In addition, the formation energy of the 6C ring is found to be 1.2 eV larger than that for the dimer (0.5 eV with the PBE[49] functional) and this is the second lowest formation energy. Considering carbon defects may be created by kinetic processes in experiments, for the binding energy discussed in Supplementary Note 1, the 6C ring has the largest binding energy[50] among the considered carbon clusters, which means the defect will agglomerate if it can diffuse. Having identified the 6C ring defect as a stable defect configuration, we now focus on its structural and electronic properties. In the neutral charge state, the ground state configuration of the defect embedded in the hBN layer is a closed-shell singlet, and it exhibits D3 symmetry. The electronic structure of hBN with the 6C ring defect is shown in Figure ; it features two pairs of degenerate e″ orbitals where two e″ orbitals fall close to the valence band maximum, fully occupied by four electrons, and the other two fall close to the conduction band minimum. The electronic configuration reads as |eo″eo″eu″eu″⟩, where o and u indicate the occupied and unoccupied states, respectively. This leads to the 1A1′ symmetry of the ground state.

Figure 1

(a) Different carbon defects we considered here and a simulated scanning tunneling microscopy image for the 6C defect. (b) Calculated formation energy vs Fermi level under N-rich and N-poor conditions. The gray color depicts the band edge.

Figure 2

(a) Single-particle energy level of the carbon ring defect in the ground state. The subscripts o and u indicate the occupied and unoccupied defect states, respectively, while the arrows denote the spin directions. (b) Wave function isosurface of defect levels. (c) Energy diagram of the optical transition with the zero-phonon line (ZPL) and Huang–Rhys (HR) factor calculated with density functional theory. The values in parentheses are the corrected ZPL with the product Jahn–Teller (pJT) effect. The right schematic figure represents the four-layer APES of the pJT effect. The dashed line is the energetically global minimum loop.

From the group theory analysis, the electronic transitions between the e orbitals give rise to four excited states in both singlet and triplet manifolds, expressed as follows:Due to the high degeneracy of the defect orbitals in D3 symmetry, each excited state represents a combination of two Slater determinants (see Supplementary Note 2). Furthermore, each of the single-electron transitions leads to the Jahn–Teller instability for both occupied and empty defect orbitals; this is achieved via a coupling to a quasi-localized E vibration mode and is known as a product Jahn–Teller (pJT) effect.[51−53] Thus, the total Hamiltonian, which accounts for both electronic correlation and pJT, is given aswhere a, a, a†, and a† are ladder operators for creating or annihilating the E phonon mode in the two-dimensional space while the first term is the vibrational potential energy of the system. Ŵ is the electronic Hamiltonian, and ĤJT is the JT part.

To solve the Ĥtot, we first construct the Ŵ. Here, the single determinants, which constitute the wave functions in eq , are shown in panels a and d of Figure . In D3 symmetry, the four single determinants form two double-degenerate branches with Ed(|eo″eu″⟩) = Ed(|eo″eu″⟩) and Ed(|eo″eu″⟩) = Ed(|eo″eu″⟩), where Ed is the total energy of the (diabatic) state. In the singlet manifold, eo″ → eu″ (or eo″ → eu″) configurations are stabilized over 41 meV by the exchange interaction [so that Ed(|eo″eu″⟩) is lower than Ed(|eo″eu″⟩)], while their order is reversed for the triplets.

Figure 3

Single-particle energy level diagram of the carbon ring defect for (a) singlet and (d) triplet excited states. The filled and empty arrows indicate the occupied and empty states with up and down spin directions, respectively. (b and e) Calculated APES for the singlet and triplet states, respectively. The dots are from DFT results, and the solid line is fitted on the basis of the pJT model. The standard deviation is <3%. X = 0 is the geometry with D3 symmetry, and the energy minima could be achieved by removing the symmetry restriction. (c and f) Energy diagrams for the four states with the TDDFT method for the singlet and triplet states, respectively. The coordinates are built on DFT optimization. The pJT effect is not included here.

Due to the complex nature of the excited states, the electronic Hamiltonian needs to be defined by using a robust method for the excited states. To this end, we compute the excitation energies of the 6C defect by CC2, focusing on a representative flake model. These calculations were assisted by TDDFT to access the transition properties, as well as by two other post-Hartree fock methods (SOS-ADC2 and NEVPT2). The resulting (vertical) excitation energies, obtained at the HSE geometry, are summarized in Table 3 of the Supporting Information. Here, we found that all of the approaches consistently predict the appearance of the localized excited states in the energy range between 4 and 5 eV. It is noteworthy that, at the high symmetry point, the two lowest A1′ and A2′ states are dark, while the transitions to E′ are optically allowed, whic is evident because of the value of the oscillator strength (∼0.93 atomic unit). From these calculations, using the definition from refs (51) and (52), the electronic Hamiltonian is expressed as follows (see Supplementary Note 3)where A1′ and A2′ are nondegenerate states and E′ is a double degenerate state. Here, Λ and Δ indicate the static electronic correlation energies for the A1,2′ and E′ states, respectively. The values obtained by TDDFT are −175.5 and −634.5 meV, respectively, for the singlets and 260.5 and 74.5 meV, respectively, for the triplets.

Having defined Ŵ, we now focus on the pJT Hamiltonian, given aswhere σ̂ = |e⟩⟨e| – |e⟩⟨e| and σ̂ = |e⟩⟨e| + |e⟩⟨e| are Pauli matrices; σ̂0 is the unit matrix, and σ̂0 = |e⟩⟨e| + |e⟩⟨e|. Fo and Fu are the electron–phonon coupling coefficients, and the major effect of the strong electron–phonon coupling is to drive the excited states out of D3 symmetry to a lower C2 symmetry by elongating two of six C–C bonds. The JT energies, denoted as EJT1 and EJT2 for |eo″eu″⟩ and |eo″eu″⟩, respectively, are determined by fitting the adiabatic potential energy surfaces (APES) from ab initio results, as shown in Figure . We found that the JT effect is much more significant for |eo″eu″⟩ than for |eo″eu″⟩, which yields the negligible EJT2. More specifically, the values of EJT1 are 187 and 239 meV for the singlets and triplets, respectively, while the EJT2 values are only 0.46 and 0.14 meV, respectively. The effective vibration energy ℏωE is then deduced from the lowest branch of the APES parabola in dimensionless generalized coordinates. The detailed construction of the pJT Hamiltonian is provided in Supplementary Note 4.

The solutions for the total Hamiltonian from eq that incorporate the vibrational and electronic parts for the singlet and triplet states are plotted in Figure . For the singlets in D3 symmetry, the states appear in the following order: E(A2′) < E(A1′) < E(E′). A2′ shows no sign of the JT instability or a mixture with E′; thus, it maintains a high-symmetry configuration and remains dark along the configuration coordinate. By contrast, when the system is driven out of D3 symmetry, the mixing between A1′ and E′ is apparent.

Figure 4

Eigenvalues for the total Hamiltonian of the system in one dimension (Y = 0) for (a) singlet and (b) triplet states. Data for pure states A1′, A2′, and E′ are depicted with black, red, and cyan dots, respectively. The lowest APES branch is a mixed state of A1′ and E′. (c) Polaronic eigenstates for the (left) singlet and (right) triplet with full rotation. The second-order pJT strength could be estimated by the energy splitting between the two lowest eigenvalues. (d) Schematic energy diagram of the electronic states and possible ISC transitions. The black dashed line links states with the same representation in different spin manifolds. The green line links states enabled by pJT-induced mixing that happens between states labeled in orange. (e) Simulated PL spectrum (red) and experimental data (black dots). The PSB of isotope 13C is also shown. The ZPL position is aligned by 0.08 eV to match the first peak in the PSB. The Gaussian broadening is 10 meV. Four peaks can be identified at 4.095, 3.905, 3.711, and 3.551 eV, which are consistent with experimental observation. The inset is the schematic coordinate diagram of the isotopic effect. (f) Simulated PL spectrum of the dimer (CNCB), 4Cpair, and 6C ring where the ZPL energies are aligned for the sake of comparison of PSBs.

A direct diagonalization of the total Hamiltonian with the pJT and electronic part ⟨Φ̃|Ĥtot|Φ̃⟩ is shown in Figure c (see Supplementary Note 5). A converged solution demonstrates that the lowest eigenstate contains 68% of the Ã1′ component in the singlet manifold (and 63% in the triplet manifold). The energy splittings between the lowest two eigenvalues are 7.1 and 3.1 meV for the singlets and triplets, respectively. On the basis of the degeneracy of polaronic levels, we assigned the lowest state to Ã1′ and the second one to Ẽ′. Given that only Ẽ′ is bright, the process requires a thermal activation and results in strongly temperature-dependent PL emission. Furthermore, the position of the ZPL based on the full Hamiltonian is calculated as followswhere Ee and Eg are the energies of excited state and ground state, respectively. The computed value is 4.21 eV, which closely agrees with the experimental data.

To further support the validity of our model calculations, we approach the A1′ geometry by TDDFT and CC2. The robust CC2 approach predicts a decrease in the symmetry to C2, while the TDDFT method preserves the D3 symmetry. Here, the energy gap between A1′ and E′ reflects the magnitude of the electronic coupling between the respective diabatic states. In the case of TDDFT, the value (223 meV) is considerably larger than that from CC2 (178 meV); this points to a strong coupling regime, where two diabats develop a single minima on the APES.[54] For the 6C defect, this relaxation is particularly important, because the coupling to the E phonon mode enables the intensity borrowing from the allowed E′; otherwise, the A1′ state remains optically forbidden.

To clarify the discrepancy between TDDFT and CC2 for the excited state geometry, we performed the TD-PBE0 calculations in a periodic monolayer. The results, shown in Table 4 of the Supporting Information, indicate that the flake model provides a reasonable description of the vertical spectrum relative to the periodic structure. However, while the difference between the A2′ energies is only 37 meV, it gradually increases to 298 meV for E′. Indeed, in this case, the electronic coupling between A1′ and E′ decreases to a much smaller value of 135 meV, and therefore, the stabilization of the C2 configuration is expected. We associate this behavior with the quantum confinement effect, which appears to be more harmful for TDDFT than for the wave function-based methods.

The simulated PL spectrum including the pJT distortion is shown in Figure e. Here, four prominent peaks in the phonon sideband with an averaged energy space of 180.3 meV perfectly match the experimental PL spectrum.[19] From these calculations, we also determine a Huang–Rhys (HR) factor, S, of 2.16, which closely agrees with the experimental results (S = 1–2). In addition, with the CC2 approach, we obtained a HR factor of 1.3 for the heteroatoms forming the flake. It is noteworthy that, at the relaxed A1′ geometry, the CC2 approach predicts that the wave function is governed by a single determinant with a relative contribution of 83%. This justifies the application of the ΔSCF for computing the vibronic sideband of A1′.

Next, we evaluate the radiative lifetimes on the basis of the following expressionwhere ϵ0 is the vacuum permittivity, ℏ is the reduced Planck constant, c is the speed of light, nD = 2.5 is the refractive index of hBN at ZPL energy EZPL, μ is the optical transition dipole moment, and η is the fraction of E′ in the polaronic state. The symmetry lowering makes Ẽ′ less bright (see Table 5 of the Supporting Information), yielding a τrad of 1.54 ns at room temperature (2 ns at 150 K for the SPE experiment[6]). This value is temperature-dependent considering the thermal occupation of Ẽ′. Nonetheless, it is very close to the observed value of ∼1.1 ns.[19]

A nonradiative transition occurs between A1′ and the lower-lying A2′ in the singlet manifold. This process could bleach the fluorescence if it is faster than the emission. In a low-temperature limit, the computed rate is 509 MHz (1.98 ns), which is slower than the radiative rate mentioned above. The optimal quantum efficiency for the defect is 52% at 300 K. However, this is influenced by the temperature, which can change the distribution between the dark and bright polaronic states, as shown in Supplementary Notes 6 and 7. We note that the nonradiative decay via phonons from the singlet A2′ toward the ground state is very slow due to the large gap between the two; thus, recombination of hot charge carriers via a two-photon absorption process is the likely process for obtaing the ground state once the electron is scattered to the dark singlet A2′ state.

Finally, after identifying the 6C defect as a promising candidate for UV emission, we compare its properties with those of 4C and CNCB. While the CNCB defect was described elsewhere,[27] for 4Cpair we computed a ZPL of ∼4.4 eV and a HR factor of 1.9. As demonstrated in Figure f, all three defects exhibit a remarkably similar phonon sideband. The minor differences between those are seen in the intensities of the replicas at the lower energies. These findings are in line with a recent experimental work in which a continuous distribution of ZPL lines around 4.1 eV[23] is observed. The similarities of PL features among these three defects indicate that other experimental techniques are needed to distinguish among those. In particular, in Supplementary Note 8, we show that the three defects demonstrate a slightly different blue shift with respect to the content of isotope 13C. Moreover, as discussed in Supplementary Note 9, we found a striking difference for these defects when considering a response to the applied strain. More specifically, the optical intensity for 6C is found to be largely affected by the uniaxial strain, but in the case of CNCB, only a weak effect is observed. Therefore, although different carbon pairs are experimentally feasible, it is essentially the 6C defect that permits a high sensitivity of the signal to the external perturbations.

In summary, on the basis of an extensive theoretical investigation, we explored the potential of substitution of carbon defects for developing UV single-photon sources in hBN. We found that carbon atoms are preferentially arranged into chains, which are stabilized to the formation of energetically favorable C–C bonds. Of those defect configurations, we identified several potential candidates for UV emission, including CNCB, 4C, and 6C defects, because they feature a photostable (neutral) charge state. The 6C defect of which configuration was observed in the experiments exhibits a highly nontrivial emission mechanism in which the second excited state is optically activated by the product Jahn–Teller effect. More specifically, the ZPL is computed at 4.21 eV and the HR factor is found to be 2.1. The simulated PL spectrum shows the phonon replicas with an energy spacing of 180 meV. The upper limit of the estimated radiative lifetime is ∼1.17 ns. All of these properties closely resemble the PL signal that is present in many hBN samples. Given the relatively low formation energy and complete agreement with the experimental measurements, these results outline the 6C defect as a plausible source of the observed UV emission. We infer that the 4.1 eV PL signal likely appears as a commutative effect from different types of point defects. Furthermore, it is likely that the 6C ring defect is responsible for the temperature[15] and strain dependency of the emission from the family of 4.1 eV emitters.