A DFT study of electron–phonon interactions for the $\text{C}_2\text{C}_\text{N}$ and $\text{V}_\text{N}\text{N}_\text{B}$ defects in hexagonal boron nitride: investigating the role of the transition dipole direction

Single-photon sources (SPSs) are essential components of many quantum technologies, including linear optical quantum computing [1, 2], quantum random number generation [3, 4], and quantum communication [5–7]. Atom-like solid-state emitters are regarded as promising SPSs because they combine the excellent optical properties of atoms with the reproducible and scalable fabrication methods that are available to solid-state materials [8–11]. Among solid-state SPSs, emitters in two-dimensional (2D) hexagonal boron nitride (h-BN) [12–14] have recently generated significant interest due to experimental observations of ultra-bright emission at room temperature [15–22].

Many quantum technologies which incorporate SPSs make use of single- or two-photon interference where indistinguishable photons are required. A significant disadvantage associated with solid-state systems such as h-BN is that phonon-induced dephasing typically results in homogeneous spectral broadening of the emission lines, which in turn, reduces the degree of photon indistinguishability [23].

When the excited electronic state of an optical transition is subject to negligible broadening mechanisms, the zero-phonon line (ZPL) is said to be Fourier transform (FT) limited, meaning the spectral width is given by the radiative decay rate [24–27]. SPSs which produce FT-limited photons are desired because they have the potential to serve as efficient sources of indistinguishable photons. Unfortunately, solid-state SPSs are typically required to operate at cryogenic temperatures to produce FT-limited emission [27–29]. The development of solid-state indistinguishable SPSs which can operate at room temperature would represent a major advancement as it would drastically reduce the difficulties associated with cryogenic cooling.

Recently, Dietrich et al [30] observed FT-limited photons from emitters in h-BN at room temperature. While the atomic structure of the emitter (or emitters) responsible is not yet known, it has been suggested that the underlying mechanism responsible for the narrow lines is an out-of-plane distortion of the emitter's orbitals and/or the transition dipole moment [31]. Support for this claim is given by analysis of the measured photoluminescence spectrum which indicates that the emitter weakly couples to in-plane phonon modes. Additionally, the ability to observe emission from a mechanically decoupled emitter requires that the emitter is on a tilted plane, suggesting that the transition dipole is perpendicular to the plane.

Motivated by the desire to develop solid-state SPSs which are capable of efficiently providing indistinguishable photons at room temperature, we have developed an approach using first-principle calculations to determine if there is a correlation between the orientation of the transition dipole moment and the magnitude of the ZPL broadening. Density functional theory (DFT) calculations have been extensively used to predict the electronic, vibrational, and optical properties of numerous solid-state emitters [32–39]. A computational methodology for determining the electron–phonon coupling associated with an optical transition has been developed by Alkauskas et al [40, 41]. The method was applied to the negatively charged nitrogen vacancy (NV⁻) center in diamond, resulting in calculated luminescence and absorption lineshapes that are in excellent agreement with experimental data. Unfortunately, the model does not predict the homogeneous broadening of the ZPL, and the width is selected to match experimental observations.

The temperature dependence of the ZPL width for the NV⁻ center has, however, been successfully described using a model that considers interactions between the excited electronic states and acoustic phonons [42]. A similar investigation was performed for the case of the silicon-vacancy center in diamond [43], where the temperature dependence of the linewidth was correctly predicted using a model of electron-phonon interactions in the excited electronic state. The key quantities required to calculate the linewidth include the phonon density of states (DOS) and the electron–phonon matrix elements. The analyses presented in [42] and [43] were based on a combination of theory and experiment. We perform a related investigation for defects in 2D h-BN without the need for experimental data, which is not yet available for the defects analyzed in this work.

To reconcile the theory presented in [42] and [43] with the experimental observations of [31], our investigation aims to determine if there is a correlation between the orientation of the transitional dipole moment and the strength of the electron-phonon interactions which are believed to ultimately dictate the width of the ZPL. As supported by the experimental observations made in [31], it seems possible that a defect with an out-of-plane transition dipole moment may be decoupled from in-plane phonon modes, and as such, the corresponding ZPL broadening may be limited to contributions from only the out-of-plane modes. The physical intuition behind our investigation is to test if the electron-phonon interactions are indeed significantly different for defects with in- and out-of-plane transition dipole moments.

We analyze the two limiting cases of the orientation of the transition dipole; when it is perpendicular to the h-BN plane and when it is parallel to the plane. The two defects studied were selected based on this criterion. Our DFT calculations reveal that the $\text{C}_2 \text{C}_\text{N}$ and $\text{V}_\text{N} \text{N}_\text{B}$ defects have transition dipoles that are in- and out-of-plane, respectively. We use density functional perturbation theory (DFPT), the goal of which is to compute derivatives of the DFT electronic energy with respect to different perturbations, to calculate and compare both the phonon DOS and the electron-phonon matrix elements of these two defects. Our findings provide no indication that an out-of-plane transition dipole moment by itself will result in the low electron-phonon coupling required to produce FT-limited lines at room temperature.

The paper is organized as follows. In section 2, we discuss the electron-phonon interactions which are responsible for spectral broadening. In section 3, we provide the details of our calculations and present the results. In section 4, we discuss the limitations associated with our analysis, along with possible extensions to our investigation. We conclude and provide an outlook in section 5.

2.1. Electron–phonon interactions

In Kohn–Sham DFT formalism [44], an interacting many-particle system is described in terms of an effective non-interacting system. The standard Hamiltonian which describes the electron-phonon coupling within DFT is obtained by expanding the Kohn–Sham effective potential in terms of nuclear displacements from their equilibrium positions, and is given by [45]

$$\begin{align}\begin{array}{*{20}{l}} {{{\hat H}_{{\mathrm{ep}}}} = \frac{1}{{\sqrt {{N_{\mathrm{p}}}} }}\sum\limits_{\begin{array}{*{20}{c}} {{\mathbf{k}},{\mathbf{q}}}\\ {mn\nu } \end{array}} {g_{mn\nu }^{(1)}} ({\mathbf{k}},{\mathbf{q}})\hat c_{m,{\mathbf{k}} + {\mathbf{q}}}^\dagger {{\hat c}_{n,{\mathbf{k}}}}({{\hat a}_{\nu ,{\mathbf{q}}}} + \hat a_{\nu , - {\mathbf{q}}}^\dagger )}\\ {{\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} + \left[ {\frac{1}{{{N_{\mathrm{p}}}}}} \right.\sum\limits_{\begin{array}{*{20}{c}} {{\mathbf{k}},{\mathbf{q}},{\mathbf{q}}\prime }\\ {mn\nu \nu \prime } \end{array}} {g_{mn\nu \nu \prime }^{(2)}} ({\mathbf{k}},{\mathbf{q}},{\mathbf{q}}\prime )\hat c_{m,{\mathbf{k}} + {\mathbf{q}} + {\mathbf{q}}\prime }^\dagger {{\hat c}_{n,{\mathbf{k}}}}}\\ {\left. {{\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} \times ({{\hat a}_{\nu ,{\mathbf{q}}}} + \hat a_{\nu , - {\mathbf{q}}}^\dagger )({{\hat a}_{\nu \prime ,{\mathbf{q}}\prime }} + \hat a_{\nu \prime , - {\mathbf{q}}\prime }^\dagger )} \right].} \end{array}\end{align} \tag{ 1 }$$

The first line of equation (1) describes, to first order in nuclear displacements, the linear coupling between an electron in band n (m) with wave vector k (k+q), and a phonon in branch ν with wave vector q (−q). The electronic (phononic) creation and annihilation operators are given by $\hat{c}_{n, \textbf{k}}^\dagger$ and $\hat{c}_{n, \textbf{k}}$ ($\hat{a}_{\nu, \textbf{q}}^\dagger$ and $\hat{a}_{\nu, \textbf{q}}$), respectively. The phonon wave vectors define a uniform grid of $N_\mathrm{p}$ points in one unit cell of the reciprocal lattice. The second and third lines of equation (1) describe the additional quadratic electron-phonon coupling that arises when expansion of the Kohn–Sham potential includes second-order nuclear displacements. The quadratic terms describe the interactions between an electron and two phonons, where one phonon is in mode ν with wave vector q and the other is in mode $\nu^{\prime}$ with wave vector $\textbf{q}^{\prime}$. The matrix elements $g_{mn\nu}^{(1)} (\textbf{k}, \textbf{q})$ and $g_{mn \nu \nu^{\prime}}^{(2)} (\textbf{k}, \textbf{q}, \textbf{q}^{\prime})$ quantify the magnitude of the linear and quadratic electron-phonon couplings, respectively.

Previous investigations [42, 46–49] into the temperature-dependent broadening of ZPL emission from solid-state defects have identified that two-phonon processes dominate at higher temperatures, but the specific transitions involved vary from defect to defect. In the negatively charged silicon-vacancy center in diamond, for example, linear interactions with phonons result in second-order elastic scattering processes that dominate the ZPL broadening at temperatures above 70 K [43]. In the NV⁻ center, on the other hand, linear interactions with phonons result in second-order inelastic Raman-type scattering processes, and quadratic interactions with phonons result in first-order elastic scattering mechanisms [42]. Both types of processes must be considered when determining the ZPL width at room temperature.

While the exact nature of the electron-phonon interactions depends on the electronic structure of the defect, all these processes share the fact that the associated interaction rates, which are calculated using Fermi's Golden Rule [50], are functions of both the phonon DOS and the electron-phonon matrix elements. We acknowledge that the current theoretical understanding of ZPL broadening is not well-developed, however, given the success of [42] and [43] at fitting theory to experimental data it is clear that the electron-phonon coupling plays a central role. Thus, calculating the electron-phonon matrix elements and the phonon DOS provides valuable insight into the processes responsible for the ZPL broadening.

The FT-limited emission reported in [30] is expected to be the result of either very small matrix elements or vanishing phonon DOS, or both. If defects with out-of-plane transition dipole moments generally result in FT-limited emission at elevated temperatures, then the corresponding electron-phonon coupling should be noticeably smaller than that of a defect center which possess an in-plane transition dipole.

To simplify the analysis, we do not consider electron–phonon interactions which result in transitions between electronic bands. We are interested in a system that displays no ZPL broadening when the temperature increases, requiring that both intraband and interband scattering processes are negligible. Additionally, phonon-assisted transitions between excited states are expected to be negligible for both defects, since the excited states are non-degenerate and well-separated in energy [34, 51]. A good starting point is therefore to determine the impact that the transition dipole direction has on only the intraband scattering rates.

The rate at which acoustic phonon elastic scattering processes occur in the excited electronic state, $\vert e\rangle$, is given by first-order Fermi's Golden Rule [42]:

$$\begin{align} \begin{split} \Gamma_\text{elastic}^\text{2-ph}(T) & = \frac{2 \pi}{\hbar} |\langle e | \hat{H}_\text{ep} | e\rangle|^2 \, \delta(E_f - E_i)\\ & = \frac{2 \pi}{\hbar} \int_0^{E_\text{cutoff}} n ( n+1 ) \Big( \,\overline{g_{ee}^{(2)}(E)} \times \rho(E)\Big)^2\mathrm dE, \end{split} \end{align} \tag{ 2 }$$

where $\hat{H}_\text{ep}$ is given by equation (1), $n = 1/(\textrm e^{E/k_\textrm B T} - 1)$ is the thermal occupation of phonons with energy E at temperature T, $k_\textrm B$ is the Boltzmann constant, $\hbar$ is the reduced Planck constant, ρ is the phonon DOS, and $\overline{g_{ee}^{(2)}(E)}$ is the average coupling strength over the phonon modes with energy E in the excited electronic state.

The elastic scattering processes involve the absorption of one phonon and the emission of another phonon with the same energy. This process will not result in a change in the electronic state but rather a transition to and from a vibrationally excited state. The Dirac delta function ensures that energy conservation is not violated. In the case of two-phonon elastic scattering, the initial and final energies (E_i and E_f ) must be the same, and hence, both phonons must have the same energy.

The integral in equation (2) is evaluated up to the maximum energy of the acoustic phonon modes (E_cutoff). The acoustic cut-off is required in [42] to properly describe the magnitude of the ZPL width. It is not yet clear why only acoustic phonon interactions should lead to ZPL broadening, and this issue remains an open question. Despite this, the model of electron interactions with acoustic phonons is able to correctly predict the temperature dependence of the ZPL width in both the NV⁻ center [42] and the silicon-vacancy center [43], and as such, it is sufficiently well-developed to serve as the theoretical basis of our analysis.

The theory of [42] states that the linewidth of the ZPL is given by the sum of an approximately temperature-independent optical decay rate (Γ₀) and the rates associated with the dominant electron-phonon scattering interactions. For the case where all phonon scattering mechanisms are negligible, except for those given by equation (2), the ZPL linewidth ($\Delta_\text{ZPL}$) is given by:

$$\begin{equation} \Delta_\text{ZPL}(T) = \Gamma_0 + \Gamma_\text{elastic}^\text{2-ph}(T) . \end{equation} \tag{ 3 }$$

We emphasize that we do not attempt to calculate the ZPL width by evaluating equation (3). For a given defect, there may be additional electron-phonon interactions that contribute to the linewidth. Instead, our focus is on analyzing a single contribution to the ZPL linewidth, and determining if it correlates with the orientation of the transition dipole moment. To that end, we examine the first-order intraband scattering rates, which depend on the phonon DOS and the quadratic electron-phonon matrix elements, as described by equation (2).

The temperature dependence of the ZPL width given by equation (3) is determined by properties of the excited electronic state of the optical transition. The transition dipole moment, on the other hand, is a function of both the ground and excited states. Therefore, the orientation of the transition dipole moment can only play an indirect role in ZPL broadening through the properties of the excited state. That said, the observations reported in [31] clearly suggest that the decoupled emitters have transition dipole moments that are perpendicular to the h-BN plane. It is, therefore, worthwhile to determine if there is a correlation between the orientation of the transition dipole and the electron–phonon coupling that ultimately leads to ZPL broadening.

2.2. Phonon DOS

In the Born–Oppenheimer approximation, the nuclei move in a potential energy (U) given by the total energy of the electronic system calculated at fixed nuclei. Dynamical properties of solid-state structures can be derived by expanding U with respect to small displacements, $u_{\kappa \alpha}^l$, about the equilibrium configuration of the nuclei, where α are the cartesian indices of nucleus κ in unit cell l. In the harmonic approximation, the expansion is given by

$$\begin{equation} U = U_0 + \frac{1}{2} \sum_{l \kappa \alpha} \sum_{l^{\prime} \kappa^{\prime} \alpha^{\prime}} \frac{\partial^2 U}{\partial u_{\kappa \alpha}^l \partial u_{\kappa^{\prime} \alpha^{\prime}}^{l^{\prime}} } u_{\kappa \alpha}^l u_{\kappa^{\prime} \alpha^{\prime}}^{l^{\prime}} . \end{equation} \tag{ 4 }$$

Here, U₀ is the total energy when the nuclei are in their equilibrium positions, and $\partial^2 U / \partial u_{\kappa \alpha}^l \partial u_{\kappa^{\prime} \alpha^{\prime}}^{l^{\prime}} = C_{\kappa \alpha, \kappa^{\prime} \alpha^{\prime}}^{l, l^{\prime}}$ is the matrix of interatomic force constants (IFCs) [52, 53]. The classical motion of the nuclei is determined by Newtonian physics, which amounts to solving the generalized eigenvalue problem [54],

$$\begin{equation} \sum_{\kappa^{\prime} \alpha^{\prime}} D_{\kappa \alpha, \kappa^{\prime} \alpha^{\prime}} (\textbf{q} ) e_{\kappa^{\prime} \alpha^{\prime}, \nu}(\textbf{q}) = \omega_{\nu, \textbf{q}}^2 e_{\kappa \alpha, \nu} (\textbf{q} ) . \end{equation} \tag{ 5 }$$

The dynamical matrix $D_{\kappa \alpha, \kappa^{\prime} \alpha^{\prime}} (\textbf{q} )$ is defined as the FT of the IFC matrix

$$\begin{equation} D_{\kappa \alpha, \kappa^{\prime} \alpha^{\prime}} (\textbf{q} ) = \frac{1}{\sqrt{m_\kappa m_{\kappa^{\prime}}}} \sum_{l^{\prime}} C_{\kappa \alpha, \kappa^{\prime} \alpha^{\prime}}^{\,l = 0, l^{\prime}} \textrm e^{\mathrm{i} \textbf{q} \cdot \boldsymbol{R}_l } , \end{equation} \tag{ 6 }$$

where m_κ denotes the mass of nuclei κ, R _l is the position vector of unit cell l, and q is the phonon wave vector. The solutions to equation (5) are the polarization vectors $e_{\kappa \alpha, \nu} (\textbf{q})$ and frequencies $\omega_{\nu, \textbf{q}}$ corresponding to the phonon mode ν [45, 55].

Calculation of the phonon frequencies given by equation (5) on a grid of q-vectors allows for the determination of the phonon DOS, $\rho(\omega)$, which is defined as the number of phonon modes per unit of frequency. The total number of modes with frequencies between ω and $\omega + \mathrm d\omega$ is given by $\rho(\omega) \mathrm d \omega$. In a system of N nuclei, there are 3N modes, hence, $\int \rho(\omega) \mathrm d\omega = 3N$. After calculating the phonon modes associated with a given system, it is then possible to determine the corresponding electron–phonon coupling.

2.3. Electron–phonon matrix elements

The linear electron–phonon matrix elements are computed in DFPT as [56]

$$\begin{equation} g_{mn\nu}^{(1)} (\textbf{k}, \textbf{q}) = \sqrt{\frac{\hbar}{2 m_0 \omega_{\nu, \textbf{q}}} } \langle \psi_{m, \textbf{k} + \textbf{q}} | \partial_{\nu \textbf{q}} V | \psi_{n, \textbf{k}} \rangle , \end{equation} \tag{ 7 }$$

which quantifies the coupling between the Kohn–Sham electronic states $\psi_{n, \textbf{k}}$ and $\psi_{m, \textbf{k} + \textbf{q}}$. Here, $\partial_{\nu, \textbf{q}} V$ is the change in the Kohn–Sham potential associated with a phonon of wave vector q, branch index ν with frequency $\omega_{\nu, \textbf{q}}$, $\hbar$ is the reduced Planck's constant, and m₀ is the mass of the nuclei in the unit cell [57].

The DFPT software we use in section 3 has not yet implemented the ability to calculate the quadratic electron-phonon matrix elements. A widely accepted approximation [42, 46, 47] is that the quadratic elements are proportional to products of the linear elements, i.e, $g_{m n \nu \nu^{\prime} }^{(2)} \propto g_{m n \nu}^{(1)} \times g_{m n \nu^{\prime}}^{(1)}$. We assume $g_{m n \nu \nu^{\prime}}^{(2)} = g_{m n \nu}^{(1)} \times g_{m n \nu^{\prime}}^{(1)}$ and discuss in section 3.3 the impact that this approximation has on our results.

To shed light on the relationship between electron–phonon coupling and the transition dipole direction associated with solid-state defects in 2D materials, we investigate two defects induced on 2D h-BN structures: $\text{V}_\text{N} \text{N}_\text{B}$ [34, 58] and $\text{C}_2 \text{C}_\text{N}$ [51] which are depicted in figure 1. The primary reason for selecting these defects is that our calculations reveal $\text{V}_\text{N} \text{N}_\text{B}$ and $\text{C}_2 \text{C}_\text{N}$ have out-of-plane and in-plane transition dipoles, respectively. Details of both the relaxation and the transition dipole calculations are provided in section 3.1.

Zoom In Zoom Out Reset image size

Additional motivation for comparing these two defects can be seen in the corresponding relaxed atomic configurations as shown in figure 1. We find that $\text{C}_2 \text{C}_\text{N}$ lies in-plane whereas the central nitrogen in $\text{V}_\text{N} \text{N}_\text{B}$ is displaced 0.6 Å out-of-plane. The total energy of the $\text{V}_\text{N} \text{N}_\text{B}$ defect with the displaced nitrogen was found to be 37 meV lower than the case where the defect nuclei were all in-plane (see section 3.1 for details of the relaxation calculations). This is consistent with previous findings where defect centers belonging to the $\text{V}_\text{N} \text{X}_\text{B}$ family have an out-of-plane displacement of the X impurity atom [58–60]. $\text{V}_\text{N} \text{X}_\text{B}$ defects possess a double-well potential for the out-of-plane displacement of the defect. The two wells correspond to displacements of the X impurity atom above and below the plane, and the in-plane position represents a local maximum.

The double-well potential has several characteristics that are consistent with the findings of [30] and [31]. The first is that the strong localization of the defect orbitals around the impurity results in displacement of both the nuclei and the orbitals, possibly leading to a reduction of the coupling strength between the defect and in-plane phonon modes. The second characteristic is that the double-well potential provides an explanation as to why the emission polarization is different for resonant excitation and far off-resonance excitation. If the reflection symmetry of the defect is removed by interactions with a nearby surface or defect, then one of the wells will be lower in energy and the impurity atom will shift farther away from the plane. This can result in two sets of transitions which, in general, have different transition dipoles and polarization. It is possible that off-resonant excitations followed by nonradiative decay to one of the excited states will result in an optical transition with polarization which differs from that of the resonant transition. A defect with a double-well potential is consistent with the observations reported in [19] and [30], both of which suggest the presence of a second excited state which can be excited indirectly, and has a polarization that differs from that of the resonant transition. Finally, a characteristic of a double-well potential is that various spectral properties, such as spectral diffusion, could depend on the excitation process. An excitation via a different excited state could be an explanation for different susceptibility to spectral diffusion depending on resonant or off-resonant excitation, which is consistent with experimental observations [30, 31]. Thus, studying the $\text{V}_\text{N} \text{N}_\text{B}$ defect allows for a simultaneous investigation of the electron-phonon coupling associated with the case of a perpendicular transition dipole moment, and the case where the central defect nucleus and the defect orbitals are displaced from the h-BN layer.

The electronic energy levels of the ground and first excited states for $\text{C}_2 \text{C}_\text{N}$ and $\text{V}_\text{N} \text{N}_\text{B}$ are illustrated in figure 2. The electronic structures of these defects have some similarities. First, there are three molecular orbitals which are located within the band gap. Second, in the ground state configuration the lowest-lying level is fully occupied, the middle level is half-occupied, and the upper level is empty. Third, upon excitation, the spin-down electron in the lower level is promoted to the middle level. For both defects, we are interested in emission corresponding to the electronic transition of the spin-down electron from the middle to the lower level. Detailed investigations into the electronic structure of $\text{V}_\text{N} \text{N}_\text{B}$ and $\text{C}_2 \text{C}_\text{N}$ can be found in [34] and [51], respectively.

Zoom In Zoom Out Reset image size

The QUANTUM ESPRESSO [61–63] software package is a well-established computational toolkit for predicting electronic and phonon properties of defects in solid-state systems. The following DFT and DFPT [56, 57, 64] calculations and post-processing were performed using QUANTUM ESPRESSO. Various assumptions and simplifications made in our analysis are specific to the current capabilities offered by this package.

3.1. Electronic structure calculations

Our DFT calculations in QUANTUM ESPRESSO were performed using the generalized gradient approximation (GGA), with the Perdew–Burke–Ernzerhof (PBE) functional [65] for the electron exchange and correlation potentials. We used a plane-wave basis set with a wave function kinetic energy cutoff of 60 Ry and projector augmented wave (PAW) pseudopotentials [66]. The total energy was found to converge at approximately 40 Ry, however, the larger cutoff was selected to ensure the accurate calculation of the phonon modes.

Our super-cells consisted of 50 atoms corresponding to $5 \times 5$ unit cells of mono-layer h-BN. The structures were relaxed until the force on each atom was less than $1 \times 10^{-6}$ Ry/au, and we used a convergence threshold for the self consistency calculations of $1\times 10^{-12}$ Ry. Relaxation was performed by allowing the in-plane lattice vectors to relax, and the out-of-plane lattice vector was held fixed to impose a vacuum spacing of 15 Å. The electronic and phonon calculations were run within the 2D framework offered by QUANTUM ESPRESSO, where the Coulomb interaction is truncated in the z-direction [67]. Other works that also inspired the settings of our DFT calculations and serve as important references for the calibration of the pristine h-BN calculations and the defect study cases are [68–75].

The total energy converged with respect to the Brillouin zone sampling density, and the lower-limit necessary to achieve convergence, a $2 \times 2 \times 1$ k-grid, was used. Both the super-cell size and the sampling density are important factors which dictate the run-time associated with the DFPT calculations. In particular, the matrix element calculations are computationally expensive in terms of run-time. The use of larger super-cells and/or more dense k-grids resulted in calculations that were infeasible. Thus, the first step was to verify that our 50-atom super-cell calculations managed to capture several key properties of the defect that are of interest in our study.

In figure 3, we plot the electronic band structure and wave functions obtained for both defects. Our calculations yield three defect states within the energy gap as expected, and inspection of the wave functions reveal that the relative ordering of the states is the same as obtained in [34] and [51].

Zoom In Zoom Out Reset image size

Point defects in crystals can be thought of as a potential well in which one or more electrons are strongly localized. If the potential is sufficiently strong such that the kinetic energy is negligible, then the corresponding energy expectation value is independent of the wave vector. Therefore, the energy of a strongly localized defect state will be constant throughout the Brillouin Zone. Consistent with out expectations, the electronic bands of the defect states shown in figure 3 are relatively flat, that is, the energy does not vary significantly along the Brillouin zone path.

If the super-cell used in the DFT calculations is not sufficiently large, then the periodic replicas of the defect, which arise due to the periodic boundary conditions utilized in the calculations, can lead to fictitious interactions between the defect and the replica defects. These interactions can lead to a de-localization of the defect electron wave function, and the resulting band structure is no longer flat throughout the Brillouin Zone [76]. We found that the use of smaller super-cells in our calculations resulted in dispersive bands across the Brillouin zone.

To further validate the use of the $5 \times 5 \times 1$ super-cell in our calculations, we compared the calculated electronic band structure to that of a $7 \times 7 \times 1$ super-cell and conclude that the differences are negligible. The results of this comparison are shown in appendix A.

Given that the main features of the electronic band structure and defect wave functions associated with the 50-atom cell are consistent with the findings of [34] and [51], we conclude that the strong localization of the defect wave functions around the defect nuclei makes it possible to provide a reasonable description of the electronic structure using a medium size super-cell.

The two defects were initially of interest to us because inspection of the molecular orbitals, following the method presented in [77], indicated that the transition dipole moment associated with the first excited state to ground state transition of $\text{C}_2 \text{C}_\text{N}$ would lie in the h-BN plane (taken to be the x–y plane), and that of $\text{V}_\text{N} \text{N}_\text{B}$ would have transition dipole perpendicular to the plane. Inspection of figure 3 reveals this to be the case since the molecular orbitals of $\text{C}_2 \text{C}_\text{N}$ are both antisymmetric with respect to in-plane mirror symmetry, resulting in a transition dipole with a z-component that is equal to zero. On the other hand, the b₂ and a₁ orbitals of $\text{V}_\text{N} \text{N}_\text{B}$ are antisymmetric and symmetric, respectively, which indicates that the z-component of the transition dipole can be nonzero.

To quantify the direction of the transition dipoles, the WFCK2R.X module in QUANTUM ESPRESSO [61–63] was used to produce the real-space wave functions, from which the transition dipole moments were directly computed according to

$$\begin{equation} \boldsymbol{\mu} = \langle \psi_f | e \boldsymbol{r} | \psi_i\rangle, \end{equation} \tag{ 8 }$$

where ψ_i and ψ_f are the initial and final wave functions, e is the elementary charge, and r is the real-space unit vector. The results of the calculations are presented in table 1.

Table 1. Calculated transition dipole moment between the ground and first excited states of the $\text{C}_2 \text{C}_\text{N}$ and $\text{V}_\text{N} \text{N}_\text{B}$ defects. The transition dipole moment of the $\text{C}_2 \text{C}_\text{N}$ defect was determined to lie predominantly in the h-BN plane, and the transition dipole of the $\text{V}_\text{N} \text{N}_\text{B}$ was determined to lie predominantly perpendicular to the h-BN plane.

Defect	µ	µ (e Å)	$\sqrt{\mu_x^2 + \mu_y^2} / \mu$	$\mu_z / \mu$
$\text{C}_2 \text{C}_\text{N}$	$\langle b_2|e \boldsymbol{r}|a_2\rangle$	0.93	1.0	$1.3 \times 10^{-8}$
$\text{V}_\text{N} \text{N}_\text{B}$	$\langle a_1|e \boldsymbol{r}|b_2\rangle$	$2.5 \times 10^{-3}$	$2.4 \times 10^{-6}$	1.0

The transition dipole associated with $\text{V}_\text{N} \text{N}_\text{B}$ is significantly smaller than that of $\text{C}_2 \text{C}_\text{N}$. The reason for this is that the $\text{C}_2 \text{C}_\text{N}$ charge density in both the ground and excited states is much more localized to the carbon atoms, as compared to the ground state charge density of $\text{V}_\text{N} \text{N}_\text{B}$ which is spread around the vacancy.

The last column in table 1 presents the ratio of the magnitude of the out-of-plane component of the transition dipole to the total magnitude. Consistent with the graphical estimation, the transition dipole associated with the ZPL transition of $\text{C}_2 \text{C}_\text{N}$ is found to lie in-plane, and that of $\text{V}_\text{N} \text{N}_\text{B}$ is perpendicular to the plane.

3.2. Phonon calculations

3.3. Electron-phonon matrix elements

The phonon DOS for both defects were calculated following the method outlined in section 3.2 and the results are shown in figure 4. Comparing the DOS of both defects, there are no features which immediately suggest that the phonon densities are suppressed in $\text{V}_\text{N} \text{N}_\text{B}$. According to the theory of [42], contributions from low-energy acoustic phonon modes have the dominant impact on ZPL broadening. The peak at approximately 40 meV is likely to be the acoustic mode cutoff energy (the highest possible energy of the acoustic modes). Both defects have non-negligible contributions to the DOS from phonons with energies below the cutoff, indicating that if the direction of the transition dipole is significantly correlated with the electron–phonon interaction rates, then such a mechanism must manifest itself in the matrix elements.

Zoom In Zoom Out Reset image size

The calculated average electron–phonon coupling strengths for the excited state intraband processes were plotted as a function of energy and the results of both defect structures are shown in figure 4. Since it remains an open question as to why only the acoustic modes should contribute to ZPL broadening, we have plotted the entire range of phonon energies. As with the case of the phonon DOS, there are obvious differences between the results of $\text{C}_2 \text{C}_\text{N}$ and $\text{V}_\text{N} \text{N}_\text{B}$, however, there is no clear reduction in the strength of the electron–phonon coupling in the low energy acoustic mode regime. In fact, the matrix elements associated with the perpendicular transition dipole are slightly larger than the case of the parallel dipole in the acoustic regime. The matrix elements were calculated at several k-points in addition to the Γ point and the corresponding results show negligible differences when compared to figure 4.

In order to exclude the theoretical possibility that the coupling might be low for both defects, we compare the matrix elements to those associated with the NV⁻ center in diamond [80]. For phonon energies less than 20 meV, the matrix elements associated with the $\text{C}_2 \text{C}_\text{N}$ and $\text{V}_\text{N} \text{N}_\text{B}$ defects are approximately an order of magnitude greater than those associated with the NV⁻ center, and up to two orders of magnitude greater at energies above 20 meV. Further, our DFT calculations indicate that the DOS associated with the NV⁻ center is less concentrated at lower energies as compared to the two defects studied here. The NV⁻ center calculation was performed using a super-cell consisting of 63 atoms, compared to the 50 atom super-cell of $\text{C}_2 \text{C}_\text{N}$ and the 49 atom super-cell of $\text{V}_\text{N} \text{N}_\text{B}$. Since there are 3N modes, the NV⁻ center has more modes, however, the phonon DOS is still less concentrated at lower energies. Since the NV⁻ center is known to have ZPL widths on the order of 10³ GHz at room temperature [42], the fact that the defects studied here have larger matrix elements and DOS at low phonon energies as compared to the NV⁻ center suggests that it would be unreasonable to expect that either of these two defects would exhibit FT-limited linewidths at room temperature.

The lack of a clear decoupling in either the DOS or the matrix elements leads to the conclusion that the direction of the transition dipole alone is not sufficient to result in an appreciable lowering of the electron-phonon coupling associated with defects in 2D materials. Similarly, the out-of-plane displacement of the central nitrogen atom in $\text{V}_\text{N} \text{N}_\text{B}$ is not expected to be, by itself, the mechanism responsible for the FT-limited lines reported in [30] and [31]. Finally, we see no indication that the orientation of the transition dipole can explain the low-energy gap in the one-phonon bands presented in [31].

It is worth mentioning that our analysis of the electron-phonon coupling is not directly related to the phonon sideband (PSB). The calculated couplings presented here correspond to intraband electronic transitions. The PSB, on the other hand, involves interband transitions, namely, transitions between the excited and ground states. An interesting future project would be to use related DFT calculations to calculate the PSB of a given defect. We expect that for defects in h-BN, we should see relatively large couplings between the excited and ground states at approximately 160 meV since this feature appears in many of the experimental observations [36].

The results presented in the preceding section are subject to a number of assumptions and limitations. We suspect that the 50-atom super-cells used to model the defects represent the largest shortcoming associated with our analysis. The electronic calculations revealed that periodic replicas of the defect are sufficiently far from one another to adequately describe the localized charge density of the defect states, however, we expect that the results of the phonon calculations would be slightly different from the case where the defect is embedded in pristine h-BN with no neighboring defects. The hypothesis that a perpendicular transition dipole will result in low electron-phonon coupling is still testable if there are defects close to the defect under investigation, which is in effect what the relatively small size of our super-cell is modelling. If the direction of the transition dipole is correlated with the ZPL width, then this effect should be observable whether or not there are nearby neighboring defects. In any case, our calculations serve as a stepping stone to enable finite-size supercell corrections to be implemented that can further improve these findings quantitatively. Examples of such finite-size correction methods that can be tested in electronic structure calculations have been vastly published, e.g. [81–87]. Similar methods based on extrapolation schemes could be used to analyze finite-size effects in the phonon modes.

It is important to note that the calculations and results presented in sections 3 and 4 consider the case of a defect embedded in a single layer of h-BN. A vacuum spacing of 15 Å in the out-of-plane direction was selected such that the periodic replicas of the h-BN layers are effectively isolated from each other, and as such, the van der Waals interaction is not included in the calculation. As shown in [58], when the interlayer van der Waals interaction is included in a calculation of the $\text{V}_\text{N} \text{N}_\text{B}$ defect embedded in bulk h-BN, the out-of-plane displacement of the central nitrogen is reduced to 0.42 Å, compared to 0.60 Å for mono-layer. As such, it is possible that in the electron-phonon coupling for $\text{V}_\text{N} \text{N}_\text{B}$ in bulk h-BN could differ from that of mono-layer. That said, if the direction of the transition dipole alone is sufficient to dictate the strength of the electron-phonon couplings that ultimately lead to ZPL broadening, then this effect would be observable in our calculations despite the differences in geometry of the $\text{V}_\text{N} \text{N}_\text{B}$ defect embedded in bulk and mono-layer h-BN.

The low sampling density of both the k- and q-points is another significant limitation in our analysis. Fourier interpolation allows for calculation of both the phonon DOS and the matrix elements on fine grids, however, it is possible that subtle details are not adequately captured with the coarse grids. Along the same lines, our calculations are not able to resolve the mode-specific contributions to either the matrix elements or the phonon DOS. It is possible that the in-plane acoustic mode is indeed decoupled from the electronic state, but this information is hidden in our analysis which considers the coupling of all modes.

The phonon mode index in the EPW output is labelled by increasing phonon energies. This becomes problematic when the phonon bands cross. It is not possible to determine which phonon mode is which in the output of EPW without a direct comparison of the phonon polarization vectors. While this would be feasible for a single path through the Brillouin zone, it becomes a monumental task to identify each mode across the $100 \times 100 \times 1$ q-grid. Several attempts were made to compare the coupling strengths of the in- and out-of-plane acoustic modes in q-point regions where the acoustic bands could be identified. Unfortunately, the results strongly depend on the q-path on which the coupling was evaluated. This is why EPW is required in the first place; the electron-phonon coupling varies significantly throughout the Brillouin zone, and therefore, the calculation is required to be evaluated on a very dense q-grid. As such, we suggest that future developments of the EPW package could include mode-resolved electron-phonon matrix elements.

Another limitation is that the electronic calculations did not utilize the Heyd–Scuseria–Ernzerhof (HSE) hybrid functional [88], which has been shown to provide an accurate estimate of the band gap of h-BN, because the DFPT calculations implemented in QUANTUM ESPRESSO do not support the use of hybrid functionals. The impact of not using the HSE functional can be seen in the underestimation of the energy spacing between defect states, and by extension, an underestimation of the ZPL energy. The matrix element calculations do not consider interband transitions, as discussed in section 2.1, and therefore, underestimating the energy spacing between electronic states is expected to have negligible impact on our results.

Additionally, the electronic and phonon calculations were performed in the ground-state nuclear configuration. A more accurate description of the excited electronic state and the corresponding phonon modes could be provided by using spin-polarized, constrained occupancy calculations with the HSE functional. We suggest that future software developments allow such calculations to be run when determining both the phonon DOS and the electron-phonon coupling. When available, the electron–phonon coupling calculated in the excited state framework can be compared to that of the ground state framework presented here.

An extension of our analysis would be to calculate the ZPL width as a function of temperature. This is an interesting route, not only to identify the defects with narrow lines at room temperature, but also to identify solid-state defects in general which have unknown atomic structures. As demonstrated by the NV⁻ center and the silicon-vacancy center in diamond, defects in the same host material can have very different electron-phonon interactions, and the ZPL widths have different temperature dependencies. Thus, the ability to calculate the ZPL broadening would allow researchers to make additional theoretical predictions which could then be compared to experimental observations.

To calculate the ZPL width as a function of temperature, one would need the ability to directly calculate the quadratic electron–phonon matrix elements, determine the acoustic mode cutoff energies, as well as the symmetry-specific mode contributions to the phonon DOS [42]. Improvements made in these areas would represent a step closer towards a complete description of the line broadening mechanisms from first-principles.

The discovery of quantum emitters in layered h-BN structures that are decoupled from in-plane phonons reveals that defects in solid-state systems have the potential to serve as indistinguishable SPSs which can operate at room temperature. The underlying mechanism responsible for the decoupling remains a mystery, however, the leading hypothesis is that an out-of-plane distortion of the emitter's orbitals leads to a reduction of the electron–phonon coupling. Here, we have tested the theory that the direction of the transition dipole moment is correlated with the rate at which electron-phonon interactions occur. Our DFT calculations show that the transition dipole of $\text{C}_2 \text{C}_\text{N}$ is in-plane, and that of $\text{V}_\text{N} \text{N}_\text{B}$ is out-of-plane. DFPT calculations were used to investigate and compare the phonon DOS and the electron–phonon coupling associated with these two defects.

Our results show no clear indication that the electron–phonon interactions are drastically different for the two defects, suggesting that an out-of-plane transition dipole moment in 2D h-BN is not sufficient to avoid the mechanisms which are expected to be responsible for temperature-dependent ZPL broadening. Our findings also indicate that the out-of-plane displacement of the nitrogen impurity atom in $\text{V}_\text{N} \text{N}_\text{B}$ does not lead to an obvious reduction in the electron-phonon scattering processes.

We have discussed the various limitations associated with our calculations in hopes that future software developments can address these issues. Our detailed description of calculating the electron-phonon coupling will be of interest in the field of solid-state quantum information processing, and adds to the ever-growing toolbox of DFT methods available to researchers.

Since all mechanically decoupled emitters were found in bulk-like h-BN, the multi-layer nature of the host material could be an essential ingredient for the decoupling mechanism [89]. The many layers in bulk-like h-BN could provide a trap, for example, for molecules localized between the layers. In such a scenario, the acoustic coupling to in-plane phonon modes could be suppressed while optical coupling would still be present in agreement with experimental indications. Investigation of this hypothesis represents an interesting future research direction.

No new data were created or analysed in this study.

The authors thank Jordan Smith for useful discussions. This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) through its Discovery Grant, Canadian Graduate Scholarships, CREATE, and Strategic Project Grant programs; by Alberta Innovates through its Graduate Student Scholarship program, and by the Alberta government through its Graduate Excellence Scholarship (AGES) program. We also acknowledge the Advanced Research Computing (ARC) IT team of the University of Calgary, and the Digital Research Alliance of Canada (alliancecan.ca/en) for computational resources.

We performed an additional set of calculations using a larger super-cell with the goal of validating the calculated electronic structure of the 50-atom $5 \times 5 \times 1$ super-cell presented in section 3.1. Our additional calculation utilized a $7 \times 7 \times 1$ super-cell for the $\text{C}_2\text{C}_\text{N}$ defect, which consists of 98 atoms, and executed the electronic calculations utilizing the same parameters as discussed in section 3.1. The electronic band structure for both the $5 \times 5 \times 1$ and $7 \times 7 \times 1$ super-cells is shown in figure 5.

Zoom In Zoom Out Reset image size

Comparing the band structures of the two calculations, we note that the relative energy spacing of the defect states, highlighted in red, is very similar. The only significant differences are that, for the $7 \times 7 \times 1$ super-cell, the top of the valence band is slightly closer in energy to the b₂ band, and the b₂ is sightly more dispersive, as compared to the $5 \times 5 \times 1$ super-cell. We believe that these differences are negligible in our analysis since we calculate the electron-phonon coupling associated with the a₂ and $b_2^{\prime}$ bands. As such, we expect that the $5 \times 5 \times 1$ super-cell is sufficiently large for the purposes of our analysis.

A convergence test was performed for the electron-phonon matrix elements with respect to the density of the k-point grid. We compare the results of $1\times 1 \times 1$, $2\times 1 \times 1$, and $2\times 2 \times 1$ grids used for both the k- and q-point sampling. Additionally, it was determined that calculations of the $\text{C}_2 \text{C}_\text{N}$ defect were feasible in terms of computational run times for a $6 \times 6 \times 1$ k-grid if the q-point sampling was restricted to a $2\times 2 \times 1$ grid. Calculations of the $\text{V}_\text{N} \text{N}_\text{B}$ defect take significantly longer due to the non-planar structure, however, demonstrating convergence with the $\text{C}_\text{2} \text{C}_\text{N}$ defect is sufficient as both super-cells are the same size.

Figure 6 shows the results of our convergence test for the $\text{C}_\text{2} \text{C}_\text{N}$ defect. The average matrix elements corresponding to the $2\times2\times1$ k-grid for phonon energies above 50 meV are very similar to those of the $6\times 6 \times 1$ grid. There are slight discrepancies between the two calculations at lower energies, however, it is clear that the $2\times 2 \times 1$ results are converging to those of the $6 \times 6 \times 1$ grid.

Zoom In Zoom Out Reset image size