Introduction

In recent developments in the field of two-dimensional materials, nanodevices utilizing graphene and hexagonal boron nitride (hBN), have undergone a unique evolution. Initially, graphene served as the active material while hBN functioned solely as a substrate or capping layer due to its excellent chemical stability1,2,3. However, the roles have now reversed, with graphene primarily acting as an electrode while hBN has emerged as the active material for applications in optoelectronics, quantum optics, and quantum information science. The intrinsic properties of hBN have become a topic of great interest, including its potential for field-enhanced molecular sensing through strong coupling to molecular vibrations4,5, as well as its ability to host room-temperature magnetic textures when interfaced with metallic ferromagnets6. Furthermore, the large band gap of hBN can accommodate the localised levels from the deep point defects, providing a platform for the solid-state spin qubits with optically addressable spin states7,8,9,10,11,12,13. These spin qubits can be applied for quantum sensing14,15,16. To date, numerous defect-related single photon emitters (SPE) in hBN have been reported, covering a wide range of wavelengths from infrared to ultraviolet8,17,18,19. These SPEs hold great potential for quantum information processing, but their precise origin should still be identified to enable the realization of these applications.

Dozens of defect structures in hBN, such as native defects, carbon and oxygen impurities, to name a few, have been proposed as sources of SPEs8,9,20,21,22,23,24,25,26. However, achieving the coherent control of the spin state has mostly been demonstrated with boron vacancy11,27,28,29, the only spin defect that was unambiguously identified in hBN7,30,31. The lack of assigned defect structures is rooted in substantial variations of photoluminescence (PL) signals, which serve as the primary means of identification, across different samples of the material. While the local strain effects are commonly attributed to this behaviour32,33, the impact of other intrinsic two-dimensional phenomena, such as twisting, sliding, and variations in layer stacking34,35,36,37,38,39,40,41,42, on electronic properties of defects remains poorly understood. Sliding, for instance, occurs across a monolayer step, leading to a gradual change between two stacking patterns40,43. The twist angle can manipulate the strength of phonon-phonon coupling44, thereby influencing the design of SPEs with desirable sharp emission lines. There is compelling evidence that the properties of hBN are closely tied to its stacking sequence. As such, interfacial ferroelectricity in hBN has been reported, with electric polarization depending on the stacking order34,35. Flat bands have been the subject of several theoretical reports, revealing a splitting of the band edge states induced by different stacking patterns41,42. In the realm of defect emission, recent studies have shown that the zero phonon line (ZPL) of ultraviolet emission from defects is stacking-dependent45,46. Moreover, the brightness of this emission can be greatly enhanced by twisting the hBN layers and further tuned by an external electric field47. These intriguing findings serve as the motivation for our study, wherein we explore the effects of stacking and sliding on deep level emission in hBN.

This paper presents comprehensive theoretical calculations for the optical properties of common ultraviolet SPEs in hBN with a specific focus on different stacking sequences. Our results demonstrate a shift of ZPL of the carbon-dimer defect from \({{{\mbox{AA}}}}^{{\prime} }\) to AB stacking, which is consistent with experimental data indicating ZPL energies of 4.09 eV for \({{{\mbox{AA}}}}^{{\prime} }\) stacking and 4.14 or 4.16 eV for Bernal (AB) stacking46. Most strikingly, the choice of stacking sequences also affects the shape of PL spectrum by altering the sideband. This phenomenon can be attributed to variations in the interlayer electrostatic Coulomb interactions. Our calculations indicate that the optical lifetime of common defects in the regular stacking patterns is primarily determined by the in-plane components of the transition dipole moment. Coupled with the inversion symmetry observed along the out-of-plane direction, this property acts as a protective factor, shielding the defects from the impact of out-of-plane electric fields. However, when situated in an inhomogeneous environment from misaligned layers of hBN, the dipolar defects demonstrate a tendency to align themselves with the direction of the field. This behavior becomes most evident in twisted bilayers, where we observe a substantial alteration in the PL signals, in agreement with experimental observations47. Our findings provide useful insights into stacking-dependent emission from defects in hBN and offer a unique strategy to enhance the brightness and quantum efficiency of SPEs.

Results and discussion

Stability and electronic properties of different stacking sequences

According to the stacking order of the successive sheets or layers, a bilayer of hBN can emerge in five different high-symmetry stacking sequences48,49, as shown in Fig. 1a. First, we compared the relative stability of these sequences by computing the total energies and checking for the appearance of imaginary modes in the phonon calculations. This analysis revealed that the AA and AB″ stacking sequences are unstable, as indicated by the presence of imaginary phonon modes at 10.4 meV and 8.3 meV, respectively, as shown in Supplementary Fig. 1. These results are consistent with a previous study that reported higher total energies for these two stacking sequences48. In fact, the \({{{\mbox{AA}}}}^{{\prime} }\) stacking pattern is commonly observed in most synthesized hBN samples, and its properties have been extensively studied in the past decades50,51,52. Interestingly, we found that the AB stacking has a lower total energy than the conventional \({{{\mbox{AA}}}}^{{\prime} }\) stacking, and its presence has been confirmed by high-resolution transmission electron microscopy (HRTEM) imaging and by combining second harmonic generation (SHG) with PL spectroscopy45,53. Therefore, we focus our attention on the three stable polytypes, namely, \({{{\mbox{AA}}}}^{{\prime} }\), AB and \({{{\mbox{AB}}}}^{{\prime} }\). Our calculations were able to reproduce the experimental band gap of approximately 6.1 eV for \({{{\mbox{AA}}}}^{{\prime} }\) and AB. Note that this value neglects a contribution of the zero-point renormalization, which arises from electron-phonon interactions54. We also observe a substantial decrease in the band gap energy to 5.3 eV for \({{{\mbox{AB}}}}^{{\prime} }\), as shown in Fig. 1b. Aligned with the vacuum level, it becomes evident that the decrease in the band gap of \({{{\mbox{AB}}}}^{{\prime} }\) is primarily attributed to the shift in the conduction band minima (CBM), while the valence band maximum (VBM) remains unchanged (see Supplementary Fig. 2).

Fig. 1: Hexagonal boron nitride polytypes.
figure 1

a Schematics view of the five high-symmetry stacking sequences in bilayer hBN. In the AA stacking, boron (nitrogen) atoms in the bottom layers are fully aligned with boron (nitrogen) atoms in the top layer. The \({{{\mbox{AA}}}}^{{\prime} }\) stacking is essentially similar to AA, but the boron atoms are aligned with nitrogen atoms. The AB is formed by rotating the layers of the \({{{\mbox{AA}}}}^{{\prime} }\) stacking by 60 degrees. The difference between the \({{{\mbox{AB}}}}^{{\prime} }\) and AB″ stacking is that either nitrogen or boron atoms appear in the center of the honeycomb from another layer. b Band alignment diagrams for the stacking sequences from (a) showing the position of the band edges relative to the vacuum level and computed for the representative slab models. The respective band gap energies in bulk are 6.05 eV, 6.06 eV and 5.30 eV.

Full size image

4.1 eV-defect emitters in different stacking sequences

In the literature, several compelling models exist to accurately describe the defects responsible for single-photon emission in the ultraviolet region of hBN. These models include a 2C (CNCB or carbon-dimer) defect55, as well as our recent findings on carbon complexes, and a Stone-Wales (SW) defect21,56. The structures and energy diagrams of these defects are shown in Fig. 2. All these defects maintain a stable neutral charge state over an energy range that exceeds the ionization threshold. The 2C defect exhibits one occupied and one empty state within the band gap, both possessing a b2 symmetry. The occupied state primarily originates from a pz orbital of CN, while the empty state arises from a pz orbital of CB. On the other hand, the 4C defect features two additional levels within the energy gap between the b2 states, resulting in the lowest optical transition from a2 to \({a}_{2}^{{\prime} }\). In contrast to the previous defects, the 6C carbon ring exhibits a D3h symmetry, giving rise to degenerate e states within the band gap. These states are labeled as occupied \({e}_{{{{\rm{o}}}}}^{{\prime\prime} }\) and unoccupied \({e}_{{{{\rm{u}}}}}^{{\prime\prime} }\). Similar to the 2C defect, the SW defect also possesses one occupied and one empty defect level. Supplementary Fig. 3 depicts the detailed defect configurations that we modeled in the three polytypes. It is worth noting that the unique nature of AB stacking permits the appearance of defects in two nonequivalent configurations, denoted as AB1 and AB2. In AB1 stacking, carbon atoms are aligned with nitrogen atoms, whereas in AB2 stacking, carbon atoms are aligned with boron atoms. Despite occasional variations in the absolute energies, the basic electronic structure of these defects remains well-preserved across different stacking sequences (see Supplementary Fig. 4).

Fig. 2: Schematics view of the defects in hBN as possible SPE sources in the ultraviolet range.
figure 2

The bottom panels show the defect levels inside the band gap of hBN, obtained from the ground state HSE calculations. The 2C defect is also known as carbon-dimer defect, whereas SW refers to Stone-Wales defect.

Full size image

We present the calculated parameters, including ZPLs, Huang-Rhys (HR) factors, and radiative lifetimes, in Table 1. The results show that compared to the previous work on 2C55, our calculations yield lower ZPLs and longer radiative lifetimes. This difference can be attributed to the varying fraction of the Fock exchange. Importantly, we account for the two-determinant nature of the excited singlet state in our ZPL calculations through a correction term, as discussed in Supplementary Note 2. Our calculations reveal that interlayer interaction greatly impacts the PL spectrum. First, for all three carbon defects, the ZPLs increase when transitioning from \({{{\mbox{AA}}}}^{{\prime} }\) to AB stacking, while the variations are found to be system-specific. In particular, 2C shows a small change in ZPL, from 4.09 eV in \({{{\mbox{AA}}}}^{{\prime} }\) stacking to 4.24 or 4.21 eV with AB stacking, which is largely reminiscent of the experimental observations of Rousseau et al46. Notably, despite the substantially smaller band gap in \({{{\mbox{AB}}}}^{{\prime} }\) stacking, the ZPLs only experience a slight shift. Furthermore, the stacking arrangement affects the relaxation energy in the excited state. For the 2C defect, the relaxation energy is 0.21 eV and 0.24 eV for the AB1 and AB2 configurations, respectively. This suggests a relatively stronger electron-phonon coupling with the AB1 pattern and highlights the influence of stacking on the optical transitions of defects. This effect is further demonstrated by the calculated HR factors, which increase from 1.80 to 2.02 when transitioning from AB1-2C to AB2-2C. By contrast, the changes in lifetime range from 1.5 to 3.3 ns, which, despite notable variations in the transition dipole moment shown in Fig. 3a, are relatively small. This is because the increase in the dipole moment is offset by the decrease in ZPL energies. Thus, the primary effect of the stacking arrangements is the modification of the PL shape, with the narrowest signals observed in \({{{\mbox{AB}}}}^{{\prime} }\) stacking. To further illustrate the concept of sideband engineering by altering the stacking sequence, Fig. 3b depicts the sidebands of \({{{\mbox{AA}}}}^{{\prime} }\) and \({{{\mbox{AB}}}}^{{\prime} }\) for the 6C defect.

Table 1 Optical properties of ultraviolet emitters in hBN
Full size table
Fig. 3: Optical properties of defects in hBN polytypes.
figure 3

a Variations of transition dipole moments calculated for the three carbon defects in the different stacking sequences of hBN. b The simulated PL spectra of 6C defects in the \({{{\mbox{AA}}}}^{{\prime} }\) and \({{{\mbox{AB}}}}^{{\prime} }\) stacking sequences.

Full size image

Tuning defect emission with sliding the layers

We proceed to investigate the effect of sliding which is feasible between the \({{{\mbox{AA}}}}^{{\prime} }\) and \({{{\mbox{AB}}}}^{{\prime} }\) configurations. Due to the rapid recovery of a high symmetry configuration (either \({{{\mbox{AA}}}}^{{\prime} }\) or \({{{\mbox{AB}}}}^{{\prime} }\)), performing a full geometry optimization of the defected structures becomes cumbersome. Therefore, our focus is on the energy difference between defect levels for the sliding from \({{{\mbox{AB}}}}^{{\prime} }\) to \({{{\mbox{AA}}}}^{{\prime} }\) configurations, as illustrated in Fig. 4. We observe a decrease in the energy difference when the geometry is in an off-high-symmetry configuration, reaching a value within 0.25 eV. This observation suggests that the ZPLs may shift during the sliding process, as well. The evolution of the energy difference is primarily driven by the changes in the band gap of pristine hBN during sliding, as plotted in Supplementary Fig. 5. Generally, the band gap decreases with non-centrosymmetric stacking sequences, reaching a minimum when the other layer lies on a bridge site. There is an explanation for these effects, which suggests that chemical (Coulomb) interactions play a crucial role in determining the relative stability of different stacking sequences, while electron correlation softens the potential energy surface49.

Fig. 4: Energy difference between defect levels upon sliding.
figure 4

The calculated energy difference between the highest occupied and the lowest unoccupied defect levels as a function of sliding. Here, the configuration 1 corresponds to the \({{{\mbox{AB}}}}^{{\prime} }\) stacking, while configuration 10 is close to the \({{{\mbox{AA}}}}^{{\prime} }\) stacking. The insets show the respective configurations of the 2C defect along the sliding path.

Full size image

Tailored defect emission in twisted structures

Having described the properties of single-photon emitters in different polytypes, we now shift our focus to the defect-related PL observed in twisted hBN bilayers. Unlike symmetric translations, the locally inhomogeneous environment resulting from twisting lowers the defects’ symmetry and induces an out-of-plane net field. The latter phenomenon, in turn, can interact with the dipole moments of both ground and excited states, consequently modifying the emission spectra (also known as the Stark effect). More precisely, we examine the influence of these fields on the PL properties of 6C, 2C, and SW defects. In Fig. 5, we compare their PL sidebands calculated in the twisted bilayers and \({{{\mbox{AA}}}}^{{\prime} }\) stacking. It becomes evident that the response of these defects to twisting differs remarkably. This effect can be understood by analyzing the computed changes in dipole moments in the ground and excited states57. As confirmed by the results in Table 2, the degree of correlation between the changes in dipole moment and the HR factors is large, given that a complete relaxation toward the net fields is impeded by the hBN lattice. Specifically, the response of the 6C defect to twisting only exhibits marginal variation, while the magnitude of changes falls within the range of the effect observed in different stacking sequences. On the other hand, striking modifications are observed for the 2C and SW defects, largely due to their substantial variations in dipole moments upon excitation. It is worth noting that, due to the change in dielectric environment, the positions of the ZPL remain essentially stable for each defect.

Fig. 5: Moiré lattice of hBN with ultraviolet emitters.
figure 5

a The structure of the 2C defect in the twisted bilayer of hBN with the twist angle of 21.79∘. b–d The PL spectra simulated for the 2C, SW, and 6C defects, respectively, in the twisted bilayer of hBN. For the sake of comparison, the respective spectra are in \({{{\mbox{AA}}}}^{{\prime} }\) bulk are also provided.

Full size image
Table 2 Changes in HR factors upon twisting and difference in dipole moments
Full size table

Among these defect configurations, the response of the 2C defect is particularly intriguing, as its sideband modifications align remarkably well with experimental measurements (see Supplementary Fig. 11). Additionally, in the relaxed excited state geometry, we observe a distortion of the defect axis out of the basal plane by approximately 9 degrees. For the sake of reference, this distortion is approximately 3 degrees in the \({{{\mbox{AA}}}}^{{\prime} }\) bilayer and zero degrees in bulk. The out-of-plane distortion may have a dual effect on the observed PL intensity. Firstly, due to its alignment with the direction of the laser beam in a typical optical setup, a permanent component of the transition dipole moment (dz) enhances both absorption and emission. Note that under this scenario the total PL lifetime is weakly affected. Secondly, further modifications of the wavefunctions can result from an additional geometrical relaxation. By discriminating the contributions of these two effects, we obtain the enhancement of the PL intensity from the dz contribution by a factor of 8.9, once again aligning with experimental measurements. In turn, the total PL intensity remains rather stable suggesting that the observed effect is primary caused by the defect reorientation toward the net field. Therefore, based on a full agreement with the experimental results, we put forward the 2C defect as the primary source of the experimental response of the 4.1 eV emission observed in the twisted bilayers of hBN. In addition, Supplementary Fig. 12 demonstrates the persistent effect in a multilayer configuration, particularly near the interface between two twisted \({{{\mbox{AA}}}}^{{\prime} }\) sequences. This strengthens the comparison between our calculations and the experimental results.

It is important to note that Su et al.47 provided a different theoretical explanation for the modification of the optical signal in the twisted bilayers. However, despite utilizing an advanced electronic structure method, the authors did not consider the reorganisation of the ions consisting of the defects upon excitation. On the other hand, our calculations reveal the dominant role of the Stark effect, which appears to significantly impact the identification of defects in hBN. Moreover, these calculations offer insights for experimentally validating the proposed mechanism. As depicted in Supplementary Fig. 11, the interlayer twist activates vibrational modes within the energy range of ~10–100 meV. Specifically, the modes at close to 20 meV are responsible for the out-of-plane distortion. We believe that these signals could be observable in resonance Raman experiments and in low-temperature PL measurements, ultimately acting as the fingerprints for the environmental modulation of intra-defect emission in hBN.

Tunable defect emission and symmetry

The disparate response of defects to stacking and twisting prominently mirrors the electronic characteristics of the host material as well as the symmetry of the involved states. As shown in Supplementary Table III, defects typically exhibit the smallest HR factors when situated within a free-standing monolayer. Supplementary Fig. 13 elucidates that the host matrix induces a polarization effect in the defective layer, akin to the impact of an external field. The charge density plots further demonstrate a nuanced induced polarization across different stacking sequences, with magnitudes diminishing in the sequence of \({{{\mbox{AA}}}}^{{\prime} }\, >\, {{{\rm{AB}}}}\, >\, {{{\mbox{AB}}}}^{{\prime} }\). This observation closely aligns with our DFT calculations, underscoring that defects in \({{{\mbox{AB}}}}^{{\prime} }\) stacking consistently manifest the smallest HR factors. Noteworthy alterations in electron density across different stacking patterns predominantly occur in the basal plane of hBN, while the out-of-plane components cancel each other owing to inherent symmetry. The act of twisting has a dual impact. Firstly, it diminishes the symmetry of defects, creating new accessible pathways for coupling with electric fields. More specifically, the external fields amend the Hamiltonian with a perturbation term of Hper = e∣E∣r, where e is an elementary charge, ∣E∣ is the magnitude of the electric field, and r reflects the direction of the field. The application of selection rules is then extended to the matrix elements, defined as e∣E∣〈ϕ∣r∣ϕ〉. Table 3 summarizes the results for the defects under consideration. Notably, the 6C defect exhibits a predominant affinity for in-plane fields, both in the stacking sequences and twisted configurations. The 2C defect demonstrates coupling with in-plane fields in stacking sequences, as well. However, it specifically couples with fields inside the symmetry plane in twisted configurations, making tilting out of the basal plane the primary relaxation pathway. This observation provides clarity on its unique response to the twist effect.

Table 3 Symmetry-related properties of the defects under consideration
Full size table

Visible defect emitters in different stacking sequences

Importantly, the Stark effect can modulate the PL spectrum of defects irrespective of the emission wavelength. As illustrative examples, we focus on two visible emitters: a negatively charged VBON20 and CB8,58 in various stacking sequences. Supplementary Fig. 14 displays the electronic structure of VBON(-), revealing an optical transition between the A1 and B1 states. Symmetry analysis indicates that both states exclusively couple to an in-plane electric field. These findings are supported by DFT calculations, revealing a ZPL value variation of ≈0.08 eV. In contrast, CB exhibits a single occupied defect state within the gap, with the optical transition occurring from the occupied defect level to the (double-degenerate) conduction band. Only the excited E state in CB is notably responsive to the in-plane field, resulting in ZPL variations of up to 0.6 eV across different stacking sequences. These calculations directly show the applicability of modulating the optical spectrum through stacking and twisting for defects with emissions ranging from visible to ultraviolet.

In summary, we propose modifying the PL response of single photon emitters in hexagonal boron nitride by altering the stacking sequences. Our calculations indicate remarkable changes in the HR factors, with variations of up to 50% observed for certain defects with regular polytypes. The dipolar defects exhibit strong coupling to the polytype, indicating a prominent role of the Stark effect. Given the general interest in SPEs with sharp emission, the \({{{\mbox{AB}}}}^{{\prime} }\) stacking is expected to produce the narrowest PL signals. By introducing twisting, the effect can be further enhanced, leading to a complete transformation of the sideband shape and effectively increased brightness. Therefore, when comparing calculated spectra to experimental data, caution should be exercised given the remarkable impact of PL spectroscopy on defect identification. Our calculations for the 2C defect align particularly well with the available experimental data, suggesting that the 2C defect is likely the primary source of the 4.1 eV-emission. In turn, its exceptional signal variations could be exploited to monitor local rearrangements of hBN caused by stain, electric fields, and other perturbations. The current technique also improves the capability of PL measurements, enabling more effective identification of defect structures at large. Future investigations should focus on comprehending the coupling between defect spin properties and environmental changes in the different arrangements of hBN layers.

Methods

Density functional theory calculations

Our spin-polarized density functional theory (DFT) calculations were performed by the Vienna ab initio simulation package (VASP) code59,60 using a plane wave basis. Projector augmented wave (PAW) potentials61,62 were used with a cutoff energy of 450 eV. A 9 × 9 two-layered supercell model was constructed to avoid the interactions of defect with its periodic images (see Supplementary Fig. 4) and to apply the Γ-point sampling scheme. The interlayer vdW interaction was described with DFT-D3 method of Grimme63 for the dispersion correction. To accurately account for the band gap energy, we modified the screened hybrid density functional of Heyd, Scuseria, and Ernzerhof (HSE06)64 with a tuned mixing parameter α = 0.32 for the Fock exchange and the exchange correlation is described with hybrid generalized gradient approximation (HGGA). The geometry optimization and calculation of electronic properties were performed with the HSE functional in consistency with our previous studies. The convergence threshold for the forces was set to 0.01 eV Å−1. ΔSCF method65 was used to calculate excited states. Since the interlayer distance of different stacking sequences does not change significantly, we fix it in our model to the value of 3.31 Å. In addition, we constructed a bilayer configuration of 252 atoms to investigate the impact of twisting on the properties of SPE. To avoid lattice mismatch between the layers66, we selected a twist angle of 21.79∘ which falls into a region of intensified PL signal for the 4.1 eV-emission47.

Simulations of PL spectra and radiative lifetime

The PL spectrum was simulated using the Franck-Condon approximation by computing the overlap between the phonon modes in the ground and excited states (see Supplementary Note 4 and 6 for details)65,67. The phonon modes were calculated with the Perdew-Burke-Ernzerhof68 (PBE) functional, which is a widely used, reliable and time-saving approximation. The radiative lifetime was computed using the following equation:

$${\Gamma }_{{{{\rm{rad}}}}}=\frac{1}{{\tau }_{{{{\rm{rad}}}}}}=\frac{{n}_{{{{\rm{D}}}}}{E}_{{{{\rm{ZPL}}}}}^{3}{\mu }^{2}}{3\pi {\epsilon }_{0}{c}^{3}{\hslash }^{4}},$$
(1)

where ϵ0 is the vacuum permittivity, ℏ is the reduced Planck constant, c is the speed of light. The refractive index nD = 2.5 for hBN was chosen at the ZPL energy EZPL of around 4.1 eV.

Time-dependent density functional theory calculations

The changes in the dipole moment upon excitation were evaluated by the time-dependent density functional theory (TDDFT) method using the ORCA code69. To this end, a cluster model of 120 atoms was described with the cc-pVDZ basis set70 and the PBE0 density functional71. The TDDFT results also confirm that the lowest excitation of the 2C and SW defects are essentially between the single pairs of orbitals, which perfectly fits into the scope of ΔSCF method. A justification of ΔSCF method for the 6C defect is provided elsewhere56.

To enhance the comparison with experimental observations, where samples typically involve multiple layers, we conducted additional convergence tests using 4- and 8-layer supercell models, as discussed in Supplementary Note 7.