Exceptionally strong coupling of defect emission in hexagonal boron nitride to stacking sequences

Van der Waals structures present a unique opportunity for tailoring material interfaces and integrating photonic functionalities. By precisely manipulating the twist angle and stacking sequences, it is possible to elegantly tune and functionalize the electronic and optical properties of layered van der Waals structures. Among these materials, two-dimensional hexagonal boron nitride (hBN) stands out for its remarkable optical properties and wide band gap, making it a promising host for solid state single photon emitters at room temperature. Previous investigations have demonstrated the observation of bright single photon emission in hBN across a wide range of wavelengths. In this study, we unveil an application of van der Waals technology in modulating their spectral shapes and brightness by carefully controlling the stacking sequences and polytypes. Our theoretical analysis reveals remarkably large variations in the Huang-Rhys factors-an indicator of the interaction between a defect and its surrounding lattice-reaching up to a factor of 3.3 for the same defect in different stackings. We provide insights into the underlying mechanism behind these variations, shedding light on the design principles necessary to achieve rational and precise control of defect emission. This work paves the way for enhancing defect identification and facilitating the engineering of highly efficient single photon sources and qubits using van der Waals materials.


I. INTRODUCTION
In recent developments in the field of two-dimensional materials, nanodevices utilizing graphene and hexagonal boron nitride 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 stability [1][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 vibrations [4,5], as well as its ability to host room-temperature magnetic textures when interfaced with metallic ferromagnets [6].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 states [7][8][9][10][11][12].These spin qubits can be applied for quantum sensing [13][14][15].To date, numerous defect-related single photon emitters (SPE) in hBN have been reported, covering a wide range of wavelengths from infrared to ultraviolet [8,[16][17][18].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, and carbon and oxygen impurities, to name a few, have been proposed as sources of SPEs [8,9,[19][20][21][22][23][24][25].However, achieving the coherent control of the spin state has mostly been demonstrated with boron vacancy [11,[26][27][28], the only spin defect that was unambiguously identified in hBN [7,29,30].The lack of assigned defect structures is rooted in substantial variations of photoluminescence 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 behaviour, the impact of other intrinsic two-dimensional phenomena, such as twisting, sliding, and variations in layer stacking [31][32][33][34][35][36][37][38][39], on electronic properties of defects remains poorly understood.Sliding, for instance, occurs across a monolayer step, leading to a gradual change between two stacking patterns [37].The twist angle can manipulate the strength of phononphonon coupling [40], 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 order [31,32].Flat bands have been the subject of several theoretical reports, revealing a splitting of the band edge states induced by different stacking patterns [38,39].In the realm of defect emission, recent studies have shown that the zero phonon line (ZPL) of ultraviolet emission from defects is stacking-dependent [41,42].Moreover, the brightness of this emission can be greatly enhanced by twisting the hBN layers and further tuned by an exter-nal electric field [43].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 2C defect from AA to AB stacking, which is consistent with experimental data indicating ZPL energies of 4.09 eV for AA stacking and 4.14 or 4.16 eV for Bernal (AB) stacking [42].Most strikingly, the choice of stacking sequences also affects the shape of photoluminescence spectrum by altering the PL 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 photoluminescence signals, in agreement with experimental observations.[43].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.

II. RESULTS
According to the stacking order of the successive sheets or layers, a bilayer of hBN can emerge in five different high-symmetry stackings [44,45], as shown in Figure 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 stackings are unstable, as indicated by the presence of imaginary phonon modes at 10.4 meV and 8.3 meV, respectively.These results are consistent with a previous study that reported higher total energies for these two stacking sequences [44].In fact, the AA stacking pattern is commonly observed in most synthesized hBN samples, and its properties have been extensively studied in the past decades [46][47][48].Interestingly, we found that the AB stacking has a lower total energy than the conventional AA stacking, and its presence has been confirmed by HR-TEM imaging and by combining second harmonic generation (SHG) with photoluminescence spectroscopy [41,49].Therefore, we focus our attention on the three stable polytypes, namely, AA , AB and AB .Our calculations were able to reproduce the exper- imental band gap of approximately 6.1 eV for AA and AB.Note that this value neglects a contribution of the zero-point renormalization, which arises from electronphonon interactions [50].We also observe a substantial decrease in the band gap energy to 5.3 eV for AB , as shown in Figure 1b.Aligned with the vacuum level, it becomes evident that the decrease in the band gap of AB is primarily attributed to the shift in the conduction band minima (CBM), while the valence band maximum (VBM) remains unchanged.
In the literature, several compelling models exist to accurately describe the defects responsible for singlephoton emission in the ultraviolet region of hBN.These models include a 2C (C N C B ) defect [51], as well as our recent findings on carbon complexes, and a Stone-Wales (SW) defect [20,52].The structures and energy diagrams of these defects are shown in Figure 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 b 2 symmetry.The occupied state primarily originates from a p z orbital of C N , while the empty state arises from a p z orbital of C B .On the other hand, the 4C defect features two additional levels within the energy gap between the b 2 states, resulting in the lowest optical transition from a 2 to a 2 .In contrast to the previous defects, the 6C carbon ring exhibits a D 3h symmetry, giving rise to degenerate e states within the band gap.These states are labeled as occupied e o and unoccupied e u .Similar to the 2C defect, the SW defect also possesses one occupied and one empty defect level.Despite occasional variations in the absolute energies, the basic electronic structure of these defects remains wellpreserved across different stackings.In Supplementary Fig. 2, we provide detailed defect configurations that we modeled in the three polytypes.It is worth noting that we considered two nonequivalent lattice sites for AB stacking, 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.
We present the calculated parameters, including ZPLs, Huang-Rhys (HR) factors, and radiative lifetimes, in Ta-ble I. Our results show that compared to the previous work on 2C [51], 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 photoluminescence spectrum.First, for all three carbon defects, the ZPLs increase when transitioning from AA 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 AA stacking to 4.24 or 4.21 eV with AB stacking in accordance with experimental observations [42].Notably, despite the substantially smaller band gap in AB 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 electronphonon 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 Figure 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 photoluminescence shape, with the narrowest signals observed in AB stacking.To further illustrate the concept of sideband engineering by altering the stacking sequence, Fig. 3b depicts the PL sidebands of AA and AB for the 6C defect.
We proceed to investigate the effect of sliding which is feasible between the AA and AB configurations.Due to the rapid recovery of a high symmetry configuration (either AA or AB ), 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 AB to AA configurations, as illustrated in Fig. 4. We observe a decrease in the energy difference when the geometry is in an off-highsymmetry 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 depicted in Supplementary Fig. 9. Generally, the band gap decreases with non-centrosymmetric stackings, 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 stackings, while electron correlation softens the potential energy  surface [45].
Having described the properties of single-photon emitters in different polytypes, we now shift our focus to the defect-related photoluminescence observed in twisted hBN bilayers.Unlike symmetric translations, the locally inhomogeneous environment resulting from twisting induces an out-of-plane net field.This 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 photoluminescence properties of 6C, 2C, and SW defects.In Fig. 5, we compare their PL sidebands calculated in the twisted bilayers and AA stacking.It becomes evident that the response of these defects to twisting differs remarkably.However, this effect can be understood by analyzing the computed changes in dipole moments in the ground and excited states [53].As confirmed by the results in Table II, 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 zero-phonon line remain essentially stable for each defect.
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. 8b).Additionally, in TABLE II.Calculated changes in the dipole moments upon excitation (Debye) and HR factors, computed for the given defects in the AA stacking and twisted bilayers.For 2C and SW, the (vertical) ∆µ was computed by the TDDFT approach.For 6C, the TDDFT results in a zero-value of ∆µ; yet we assume a finite ∆µ due to the missing contribution of the product Jahn-Teller effect [52]. 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 AA bilayer and 0 degrees in bulk.The out-of-plane distortion may have a dual effect on the observed photoluminescence 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 (d z ) enhances both absorption and emission.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 d z 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.

Defect
It is important to note that Su et al. [43] 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 new insights for experimentally validating the proposed mechanism.As depicted in Supplementary Fig. 8a, 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 outof-plane distortion.We believe that these new signals could be observable in resonance Raman measurements and would further confirm the environmental modulation of intra-defect emission in hBN.

III. SUMMARY AND CONCLUSION
In summary, we propose modifying the photoluminescence 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 AB 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 photoluminescence 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.

IV. METHODS
Our density functional theory (DFT) calculations were performed by the Vienna ab initio simulation package (VASP) code [54,55] using a plane wave basis.Projector augmented wave (PAW) potentials [56,57] 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 and to apply the Γpoint sampling scheme.The interlayer vdW interaction was described with DFT-D3 method of Grimme [58] for the dispersion correction.To accurately account for the band gap energy, we modified the screened hybrid density functional of Heyd, Scuseria, and Ernzerhof (HSE) [59] with a mixing parameter α = 0.32 for the Fock exchange.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/ Å. ∆SCF method [60] was used to calculate excited states.Since the interlayer distance of different stackings 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 layers [61], we selected a twist angle of 21.79 • which falls into a region of intensified PL signal for the 4.1 eV-emission [43].
The PL spectrum was simulated using the Franck-Condon approximation by computing the overlap between the phonon modes in the ground and excited states.[60,62].The phonon modes were calculated with the Perdew-Burke-Ernzerhof [63] (PBE) functional, which is a widely used, reliable and time-saving approximation.The radiative lifetime was computed using the following equation: where 0 is the vacuum permittivity, is the reduced Planck constant, c is the speed of light.The refractive index n D = 2.5 for hBN was chosen at the ZPL energy E ZPL of around 4.1 eV.
The changes in the dipole moment upon excitation were evaluated by the time-dependent density functional theory (TDDFT) method using the ORCA code [64].To this end, a cluster model of 120 atoms was described with the cc-pVDZ basis set [65] and the PBE0 density functional [66].The TDDFT results also confirm that the lowest excitations 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 elsewhere [52].

FIG. 1 .
FIG. 1.(a) Schematics view of the five high-symmetry stackings 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 AA 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 AA stacking by 60 degrees.The difference between the AB 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

FIG. 2 .
FIG. 2. Schematics view of the three carbon clusters and a SW defect in hBN which we consider as possible SPE sources in the ultraviolet range.The bottom panels show the defect levels inside the band gap of hBN, obtained from the ground state HSE calculations.

FIG. 4 .
FIG.4.The calculated energy difference between the highest occupied and the lowest unoccupied defect levels of 2C and 6C as a function of sliding.Here, the configuration 1 corresponds to the AB stacking, while configuration 10 denotes the AA stacking.The insets show the respective configurations of the 2C defect along the sliding path.

FIG. 5 .
FIG. 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 AA bulk are also provided.

TABLE I .
Calculated ZPL energies after including the correction term (eV), HR factors, transition dipole moments (Debye) and radiative lifetimes (ns) for the carbon impurities and SW defect in the different stackings of hBN.