First-principles theory of the luminescence lineshape for the triplet transition in diamond NV centres

The excited state $^{3}E$ is an orbital doublet that forms an $E\otimes e$ Jahn–Teller system via coupling to e phonon modes [17, 30, 31]. The Jahn–Teller effect is dynamical, since the energy splitting between the vibronic sub-levels is larger than the barrier in the adiabatic potential energy surface $\delta \approx 10$ meV [31]. The presence of this effect is manifest in the broadening of the ZPL that follows a $\sim {{T}^{5}}$ rather than the usual $\sim {{T}^{7}}$ dependence at low temperatures [30, 31]. However, the effect is weak [31], and we will neglect it when calculating the phonon sideband. One way to judge the validity of this approximation is via an a posteriori comparison [17]. If a linear model of electron–phonon interactions, such as the one employed in this work, accurately describes the lineshape of the optical transition, then the Jahn–Teller effect can be considered negligible for this particular transition. As we show below, this turns out to be the case for the triplet luminescence at NV centres.

We also assume that the transition dipole moment ${{\vec{\mu }}_{eg}}$ between the excited and the ground state depends weakly on lattice parameters (the Franck–Condon approximation). At T = 0 K the absolute luminescence intensity $I\left( \hbar \omega \right)$ (i.e., photons per unit time per unit energy) for a given photon energy $\hbar \omega$ and for one emitting centre is given by (in SI units) [32]:

$$\begin{eqnarray}&&I\left( \hbar \omega \right)=\frac{{{n}_{D}}{{\omega }^{3}}}{3{{\varepsilon }_{0}}\pi {{c}^{3}}\hbar }{{\left| {{{\vec{\mu }}}_{eg}} \right|}^{2}}\mathop{\sum }\limits_{m}{{\left| \left\langle {{\chi }_{gm}} \right|{{\chi }_{e0}}\rangle \right|}^{2}}\delta \left( {{E}_{ZPL}}-{{E}_{gm}}-\hbar \omega \right).\end{eqnarray} \tag{ 1 }$$

Here ${{n}_{D}}=2.4$ is the refractive index of diamond; ${{\chi }_{e0}}$ and ${{\chi }_{gm}}$ are vibrational levels of the excited and the ground state; ${{E}_{gm}}$ is the energy of the state ${{\chi }_{gm}}$, being the sum over all vibrational modes k, i.e., ${{E}_{gm}}={{\sum }_{k}}{{n}_{k}}\hbar {{\omega }_{k}}$; and ${{n}_{k}}$ is the number of phonons of type k in this state. The absolute angle-averaged value of ${{\mu }_{eg}}$ is $\sim 5.2$ Debye, as extracted from the radiative lifetime $\tau =13$ ns of the ${{m}_{s}}=0$ spin state of the $^{3}E$ manifold [1]. A prefactor ${{\omega }^{3}}$ in equation (1) arises from the density of states of photons that cause the spontaneous emission ($\sim {{\omega }^{2}}$), and the perturbing electric field of those photons ($|\vec{\mathcal{E}}{{|}^{2}}\sim \omega$). This prefactor has to be taken into account when determining parameters pertaining to the luminescence lineshape, and this will be discussed in section 2.2. Since in both the excited and the ground electronic state the system has ${{C}_{3v}}$ symmetry, only fully symmetric ${{a}_{1}}$ phonons contribute to the sum in equation (1).

The experimental determination of the absolute luminescence intensity given in equation (1) is difficult. Thus, in this work we will consider the normalized luminescence intensity, defined as

$$\begin{eqnarray}&&L\left( \hbar \omega \right)=C{{\omega }^{3}}A\left( \hbar \omega \right),\end{eqnarray} \tag{ 2 }$$

where

$$\begin{eqnarray}&&A\left( \hbar \omega \right)=\mathop{\sum }\limits_{m}\mid \langle {{\chi }_{gm}}\mid {{\chi }_{e0}}\rangle {{\mid }^{2}}\delta \left( {{E}_{ZPL}}-{{E}_{gm}}-\hbar \omega \right),\end{eqnarray} \tag{ 3 }$$

is the optical spectral function, and C is the normalization constant: ${{C}^{-1}}=\int A\left( \hbar \omega \right){{\omega }^{3}}{\rm d}\left( \hbar \omega \right)$. $I\left( \hbar \omega \right)$ is related to $L\left( \hbar \omega \right)$ via $I\left( \hbar \omega \right)={{n}_{D}}/\left( 3C{{\varepsilon }_{0}}\pi {{c}^{3}}\hbar \right)L\left( \hbar \omega \right)$.

The evaluation of the overlap integrals $\left\langle {{\chi }_{gm}} \left| \right. {{\chi }_{e0}}\right\rangle$ immediately poses a challenge. Vibrational modes that enter into equation (3) are not those of the pristine bulk, but rather those of the solid with a defect. The use of bulk modes can lead to large discrepancies with experiment, as we will show in section 4. Lattice imperfections induce localized or quasi-localized vibrational modes that depend on the local electronic structure; in addition, the normal modes in the excited state and the ground state can be in principle quite different [33]. This results in highly multidimensional integrals that can in practice be evaluated only for molecules [34], small atomic clusters [35], or model defect systems [36].

Some kind of approximation is thus unavoidable. Here we assume that (i) the normal modes that contribute to the luminescence lineshape are still those of the solid with a defect, but (ii) the modes in the excited electronic state are identical to those in the ground state. Such an assumption is implicit in virtually all studies of defects in solids [37, 38]. First-principles calculations [19, 20], as well as a comparison of experimental absorption and emission spectra [16], indicate that the assumption does not strictly hold for the NV centre. Since the more exact calculation is not feasible, the validity of this approximation has to be checked by comparing the results with the experimental spectrum.

When vibrational modes in the ground and the excited state are identical, the optical spectral function $A\left( \hbar \omega \right)$ (equation (3)) can be calculated using a generating function approach proposed by Lax [38], as well as Kubo and Toyozawa [39]. The fundamental quantity that has to be calculated is the spectral function (also called spectral density) of electron–phonon coupling [40]

$$\begin{eqnarray}&&S\left( \hbar \omega \right)=\mathop{\sum }\limits_{k}{{S}_{k}}\delta \left( \hbar \omega -\hbar {{\omega }_{k}} \right),\end{eqnarray} \tag{ 4 }$$

where the sum is over all phonon modes k with frequencies ${{\omega }_{k}}$, and ${{S}_{k}}$ is the (partial) HR factor for the mode k. It is defined as [37]

$$\begin{eqnarray}&&{{S}_{k}}={{\omega }_{k}}q_{k}^{2}/\left( 2\hbar \right),\end{eqnarray} \tag{ 5 }$$

with

$$\begin{eqnarray}&&{{q}_{k}}=\mathop{\sum }\limits_{\alpha i}m_{\alpha }^{1/2}\left( {{R}_{e;\alpha i}}-{{R}_{g;\alpha i}} \right)\Delta {{r}_{k;\alpha i}}.\end{eqnarray} \tag{ 6 }$$

α labels atoms, $i=\left\{ x,y,z \right\}$, ${{m}_{\alpha }}$ is the mass of atom α (carbon or nitrogen, average atomic masses of naturally occurring isotopes were used), ${{R}_{\left\{ e,g \right\};\alpha i}}$ is the equilibrium position in the initial (excited) and the final (ground) excited state, and $\Delta {{r}_{k;\alpha i}}$ is a normalized vector that describes the displacement of the atom α along the direction i in the phonon mode k. One can use an alternative expression for ${{q}_{k}}$:

$$\begin{eqnarray}&&{{q}_{k}}=\frac{1}{\omega _{k}^{2}}\mathop{\sum }\limits_{\alpha i}\frac{1}{m_{\alpha }^{1/2}}\left( {{F}_{e;\alpha i}}-{{F}_{g;\alpha i}} \right)\Delta {{r}_{k;\alpha i}},\end{eqnarray} \tag{ 7 }$$

where ${{F}_{e;\alpha i}}-{{F}_{g;\alpha i}}$ is the change of the force on the atom α along the direction i for a fixed position of all atoms when the electronic state of the defect changes from $^{3}E$ to $^{3}{{A}_{2}}$ . The latter equation directly follows from the relationship $\left( {{{\vec{R}}}_{e}}-{{{\vec{R}}}_{g}} \right)=-{{\hat{H}}^{-1}}\left( {{{\vec{F}}}_{e}}-{{{\vec{F}}}_{g}} \right)$, where $\hat{H}$ is the Hessian matrix, different from the dynamical matrix only because of additional mass prefactors in the latter. The two formulations are completely equivalent in the harmonic approximation. In appendix A we show that if the dynamical Jahn–Teller effect is neglected the anharmonicities are indeed minute. While being in principle equivalent, the use of equation (7) instead of equation (6) offers a huge advantage when dealing with large systems, i.e., when extrapolating $S\left( \hbar \omega \right)$ to the dilute limit, and this is discussed in section 2.3 and appendix B.

Once $S\left( \hbar \omega \right)$ is determined, the spectral function $A\left( \hbar \omega \right)$ (equation (3)) is given as the Fourier transform of the generating function G(t) [38, 39]:

$$\begin{eqnarray}&&A\left( {{E}_{ZPL}}-\hbar \omega \right)=\frac{1}{2\pi }\int _{-\infty }^{\infty }G\left( t \right){{e}^{i\omega t-\gamma |t|}}{\rm d}t.\end{eqnarray} \tag{ 8 }$$

The generating function G(t) itself is defined as

$$\begin{eqnarray}&&G\left( t \right)={{e}^{S\left( t \right)-S\left( 0 \right)}},\end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray}&&S\left( t \right)=\int _{0}^{\infty }S\left( \hbar \omega \right){{e}^{-i\omega t}}{\rm d}\left( \hbar \omega \right),\end{eqnarray} \tag{ 10 }$$

and

$$\begin{eqnarray}&&S\equiv S\left( t=0 \right)=\int _{0}^{\infty }S\left( \hbar \omega \right){\rm d}\left( \hbar \omega \right)=\mathop{\sum }\limits_{k}{{S}_{k}},\end{eqnarray} \tag{ 11 }$$

is the total HR factor for a given optical transition. In equation (8) the parameter γ represents the broadening of the ZPL. In real situations this broadening has two contributions: the homogeneous broadening due to anharmonic phonon interactions [41, 42] and the inhomogeneous broadening due to ensemble averaging. Since neither of these two effects is modeled in our approach, γ is chosen to reproduce the experimental width of the ZPL.

The partial HR factor ${{S}_{k}}$ defined in equation (5) is the average number of phonons of type k emitted during an optical transition [18]. The total HR factor, defined in equation (11), is then the number of phonons off all kinds that are emitted during the same transition. The HR factor is thus an important parameter that characterizes the vibrational structure of the luminescence band. If (i) there was no additional prefactor $\sim {{\omega }^{3}}$ in the expression for the luminescence intensity in equation (2) and (ii) the vibrational modes in the excited and the ground state were indeed identical, then the weight of the ZPL would be given by [17, 37, 38, 43] ${{w}_{ZPL}}={{e}^{-S}}.$ Since this line corresponds, by definition, to zero absorbed or emitted phonons, ${{w}_{ZPL}}$ is often called the Debye–Waller factor, in analogy with x-ray scattering, where it represents the ratio of the elastic to the total scattering cross-section. Therefore, we also use this nomenclature to comply with the accepted practice.

The Debye–Waller factor ${{w}_{ZPL}}$ is a quantity that is directly measurable in experiment, and its determination is therefore unambiguous. In practical situations the HR factor is often deduced from the relationship $\tilde{S}=-{\rm ln} \left( {{w}_{ZPL}} \right)$, where we have added a 'tilde' to distinguish this quantity from the actual HR factor S, which is defined by equation (11). $\tilde{S}$ differs from S because of the additional assumption (i).

The spectral weight in $L\left( \hbar \omega \right)$ (equation (2)) moves to slightly higher energies in comparison to $A\left( \hbar \omega \right)$ due to the prefactor ${{\omega }^{3}}$. This increases the weight of the ZPL, ${{w}_{ZPL}}$, if determined from $L\left( \hbar \omega \right)$, and thus decreases the value of $\tilde{S}$ with respect to S. This distinction has to be borne in mind when comparing different experimental papers. In this paper we will consistently use ${{w}_{ZPL}}$ and S in their original definitions.

In this work the spectral function of electron–phonon coupling $S\left( \hbar \omega \right)$ (equations (4)–(6)) is calculated within DFT. The electronic, atomic, and vibrational properties of the NV centre are calculated in the supercell approach [44], whereby one defect is embedded in a sufficiently large piece of host material, which is periodically repeated. We take a conventional cubic cell with eight carbon atoms as the building block for larger supercells. The cubic supercell $N\times N\times N,$ for example, contains M = $8{{N}^{3}}$ atomic sites.

To study the electronic and vibrational structure of defects we use two different exchange-correlation (XC) functionals: the generalized gradient approximation in the form proposed by Perdew, Burke, and Ernzerhof (PBE) [45] and the screened hybrid functional of Heyd, Scuseria, and Ernzerhof (HSE) [27]. PBE is known to describe structural properties of many materials with high accuracy, but the calculated band gaps of semiconductors and insulators agree poorly with experiment, and this also affects the position of defect levels within the band gap. The HSE functional overcomes this problem by incorporating a fraction $a=1/4$ of screened Fock exchange (screening parameter $\omega =0.2$ Å⁻¹). HSE calculations yield excellent results for excitation energies for the spin-triplet optical transition in NV centres [46].

The properties of the excited state $^{3}E$ have been calculated using the constrained occupation method of Slater [47], as first applied to the NV centre by Gali et al [46]. In this method one electron from the ${{a}_{1}}$ orbital is promoted to one of the e orbitals. The electronic and the atomic structure is optimized with a hole in the ${{a}_{1}}$ state. To circumvent the problems with the Jahn–Teller distortion in the excited state, resulting from the degeneracy of the nominal ${{a}_{1}}e_{x}^{2}e_{y}^{1}$ and ${{a}_{1}}e_{x}^{1}e_{y}^{2}$ configurations, the coordinate dependence of the total energy in the excited state is studied here by constraining the configuration to ${{a}_{1}}e_{x}^{1.5}e_{y}^{1.5}$. This is a practical solution to restrict the excited state density to be the average of the two degenerate configurations, retaining a ${{C}_{3v}}$ symmetry.

We find that while the integrated parameters, for example the total HR factor S (equation (11)) converge quickly when the size of the defect supercell is increased, the convergence of the spectral function $S\left( \hbar \omega \right)$ (equation (4)) is significantly slower. This is a particular concern for the spectral function at lower energies, i.e., coupling to long-range acoustic phonons. As an example, let us consider a $2\times 2\times 2$ simple cubic supercell, containing 64 lattice sites. Without a defect, the lowest energy $\Gamma$-point vibration of such a supercell corresponds to the bulk transverse acoustic (TA) mode at the $\Lambda$ point with an energy of about 68 meV. This is even higher than the experimentally determined energy of the dominant phonon mode at the NV centre, being $\hbar {{\omega }_{0}}=63-65\;{\rm meV}.$ Clearly, supercells of this size are insufficient to determine $S\left( \hbar \omega \right)$.

To obtain converged results and determine the nature of vibrational states, we have performed calculations for a series of supercells: from $2\times 2\times 2$ (64 sites) up to $11\times 11\times 11$ (10648 sites). Since a direct approach for supercells containing more than a few hundred atoms is computationally too demanding, we have developed a special methodology to achieve this goal. First, in appendix A we show that it is an excellent approximation to calculate vibrational properties at the PBE level, since the relevant vibrational modes are very similar in the PBE as compared to the HSE functional. This presents huge computational savings, since HSE calculations are up to two orders of magnitude more expensive. Then in appendix B we present a methodology to calculate vibrational spectra and spectral functions $S\left( \hbar \omega \right)$ for very large systems. In short, the procedure is as follows. Partial HR factors for large systems are determined from equation (7). Forces ${{\vec{F}}_{\left\{ e,g \right\}}}$ in the large supercell $N\times N\times N$ ($N>3$), needed in that approach, are obtained from the calculation of a smaller supercell $(4\times 4\times 4)$ via a suitable embedding procedure, explained in appendix B. For these large defect supercells, vibrational modes and frequencies, that also appear in expression (7), have been determined by diagonalizing the dynamical matrix constructed from dynamical matrices of bulk diamond and NV centre in the $3\times 3\times 3$ supercell. The validity of the procedure relies on the fact that the dynamical matrix of diamond is rather short-ranged. Specific parameters of the procedure are determined from accurate convergence tests, and are discussed in appendix B.

Defect calculations have been performed with the vasp code [48, 49], and the interaction with ionic cores was described via the projector-augmented wave formalism [50]. A kinetic energy cutoff of 400 eV (29.4 Ry) has been used for the expansion of electronic wavefunctions. For the $2\times 2\times 2$ supercell the Brillouin zone was sampled using a $2\times 2\times 2$ k-point mesh, and $\Gamma$-point sampling was used for larger supercells.

To produce additional insights, we have also calculated $S\left( \hbar \omega \right)$ (equation (5)) and $L\left( \hbar \omega \right)$ (equation (2)) with an additional assumption, namely that phonon modes that contribute to the luminescence lineshape are those of the unperturbed host [51, 52]. For this purpose we have determined the vibrational modes of bulk diamond using density functional perturbation theory [53], reproducing earlier calculations [54]. Vibrational modes were determined on a very fine $27\times 27\times 27$ k-point grid close to the Brillouin zone centre, and a courser $9\times 9\times 9$ grid elsewhere. These calculations have been performed using the quantum espresso code [55] within the local density approximation [56]; this XC functional describes phonons modes of bulk diamond very well [54]. To evaluate $S\left( \hbar \omega \right)$, the modes were mapped to the $\Gamma$-point of the desired supercell. The contributions from the vacancy site are set to zero in equations (6) and (7), while the mass of the nitrogen atom was set to be equal to that of the carbon atom in this case.

For the $4\times 4\times 4$ supercell, the largest system for which we have performed actual electronic structure calculations, ${{E}_{ZPL}}$ was calculated to be 1.757 eV using the PBE functional, and 2.035 eV using the HSE functional. The latter is thus much closer to the experimental value of 1.945 eV. Our calculations agree with those of Gali et al [46] and Weber et al [21]. The HSE functional is clearly superior for describing the local electronic structure of the NV centre [46]. The difference of about 0.1 eV between the experimental and calculated ZPL is within the error bar of the HSE calculations, but would complicate direct comparisons between theoretical and experimental lineshapes. To enable a more meaningful comparison, in all subsequent analysis we set ${{E}_{ZPL}}$ to the experimental value. Thus, the broadening of the ZPL γ in equation (8) and the value ${{E}_{ZPL}}$ are the sole instances where information from experiment has been used in the theoretical results.

We first analyze the convergence of $S(\hbar \omega )$ when the size of the supercell is increased. In addition to providing justification for the computational procedure, such a study provides insights into the origin of vibrational modes that contribute to the phonon sideband.

In figure 2 we show $S\left( \hbar \omega \right)$ (equation (4)) and partial HR factors (equation (5)) as a function of the supercell size, from $2\times 2\times 2$ to $11\times 11\times 11$ (results for five intermediate supercells are omitted). The range of the left vertical axes for $S\left( \hbar \omega \right)$ was kept identical for all supercells, but note that this is not the case for the right vertical axes that apply to ${{S}_{k}}$. For the calculation of $S\left( \hbar \omega \right)$ δ-functions in equation (4) were replaced by Gaussians with widths $\sigma =6$ meV. HSE results are discussed here.

Zoom In Zoom Out Reset image size

In the case of the smallest $2\times 2\times 2$ supercell, only a few phonon modes contribute to $S\left( \hbar \omega \right)$. The most important of these are modes with energies 60.4 meV and 77.8 meV, in complete agreement with the results of Gali et al [19] (59.7 and 77.0 meV), who studied local vibrational modes for this size of supercell. While the energy of the first mode is close to the energy of the most pronounced phonon mode seen in experiment, this agreement is largely fortuitous, since, as mentioned in section 2.3, the lowest-energy bulk TA phonon mode in this supercell has a similar energy. The low-energy tail that represents the coupling to long-range phonons ($<$ 60 meV) is completely missing for this supercell, and the total HR factor S = 3.02 (inset of figure 2) is 20% smaller than the converged value S = 3.67.

In the case of a larger $3\times 3\times 3$ supercell, the most dominant vibration is the 45.9 meV mode. This result is an artifact resulting from the use of a small cell, since this vibration corresponds to the lowest-energy $\Gamma$-point bulk TA phonon—it is not an actual defect-derived mode. While the total HR factor increases to 3.27, $S\left( \hbar \omega \right)$ is still far from converged. This emphasizes possible dangers in drawing conclusions about local vibrational modes from small-size supercells [20].

When increasing the size of the supercell further, $S\left( \hbar \omega \right)$ slowly attains its converged form. figure 3 shows that the spectral function is essentially converged for the two largest supercells we use, even though there are still apparent changes in individual partial HR factors ${{S}_{k}}$. The peak of $S\left( \hbar \omega \right)$ occurs at $\hbar {{\omega }_{0}}=65$ meV, in excellent agreement with experimental findings (see section 4.3). This is the first time that theoretical calculations yield the energy of the peak decisively. Interestingly, the total HR factor, i.e. the integral of $S\left( \hbar \omega \right)$, is within ∼1% of the converged value already for the $4\times 4\times 4$ supercell.

Zoom In Zoom Out Reset image size

Figure 2 also allows us to draw the following conclusions about lattice distortions or, equivalently, coupling to phonons, that occur during the $^{3}E{{\to }^{3}}{{A}_{2}}$ optical transition:

• (i)  
  The 65 meV vibration is not a localized phonon mode, but a defect-induced vibrational resonance: it occurs within the spectrum of bulk phonon modes (0–167 meV). In figure 2 this result is evident from the fact that for larger supercells this mode splits into many closely spaced modes, with a simultaneous decrease of their absolute contributions. The 65 meV resonance is induced by the NV centre itself, and cannot be understood solely by considering bulk phonon spectrum. This is demonstrated in figure 3 and discussed in more detail in section 4.3.
• (ii)  
  In agreement with a general theory of vibrational broadening of luminescence lines [41], the spectral function is linear for small energies, i.e., $S\left( \hbar \omega \right)=\alpha \hbar \omega$ for $\omega \to 0$. Indeed, partial HR factors corresponding to acoustic modes scale like $1/\omega$, which, multiplied with the density of states of acoustic modes $\sim {{\omega }^{2}}$, yields this linear dependence. While this general behaviour is known [41], we emphasize that the prefactor to the linear dependence is system dependent, and only accurate atomistic calculations such as the ones presented here can provide the actual value. In our case we obtain $\alpha \approx 3.6\times {{10}^{-4}}$ me V⁻² = 360 e V⁻². Interaction via acoustic phonons has been recently proposed as a promising mechanism to couple two NV centres in nanodiamonds [57]. The coupling of isolated qubits is essential for any quantum computing protocol. Our calculations provide information about the coupling of NV centres to acoustic phonons in bulk diamond, and can be useful pursuing the ideas proposed in [57] ideas further.
• (iii)  
  99% of the lattice distortions due to the optical transition, as quantified by their contribution to $S\left( \hbar \omega \right)$, occur within $\sim 12$ $\overset{\circ}{A}$ of the NV centre. This follows from our finding that the total HR factor for the $4\times 4\times 4$ supercell is within 1% of the converged value. However, long-range relaxations, while contributing little to the total HR factor S, are manifest in the low-frequency part of $S\left( \hbar \omega \right)$, and are actually observed in the luminescence lineshape (see section 4.3).

In figure 3 we show a comparison of $S\left( \hbar \omega \right)$ calculated using three different approaches. From here on we use the following notation when we refer to our calculations: (i) 'HSE' refers to calculations where atomic displacements or forces in equations (6) and (7) are calculated using the HSE hybrid functional, but vibrational modes are calculated using the PBE functional. As discussed in section 2.3 and appendix A, calculations for smaller supercells show that vibrational modes calculated at the PBE level are very similar to HSE results. (ii) 'PBE' refers to calculations in which all quantities are determined at the PBE level. In both (i) and (ii), vibrational modes correspond to the defect system. (iii) 'Bulk phonons', refers to calculations in which atomic distortions or forces were determined at the HSE level, as in (i), but vibrational modes correspond to those of the unperturbed host. The comparison of (i) and (iii) should inform us whether the introduction of the defect modifies the vibrational spectrum, and whether the phonon sideband can be understood by considering bulk modes alone.

$S\left( \hbar \omega \right)$, calculated at the PBE level, is qualitatively very similar to the HSE result. The function has a peak at $\hbar \omega =64$ meV, but the absolute value of $S\left( \hbar \omega \right)$ is smaller for almost all energies. In particular, the total HR factor is $S=2.78,$ a quarter smaller than in HSE. In contrast, when the bulk phonon spectrum is used, $S\left( \hbar \omega \right)$ is even qualitatively completely different. In this case the spectral function closely follows the density of vibrational states of bulk diamond [54, 58], with a pronounced peak at $\hbar \omega \;\approx$ 150 meV. The total HR factor is 4.48 in this case. However, the coupling to low-energy ($<$ 20 meV) acoustic modes is very similar to cases (i) and (ii); indeed long-range phonons are expected to be little affected by the presence of the defect.

In figure 4 we compare the luminescence lineshape $L\left( \hbar \omega \right)$ (equations (2) and (8)), calculated using the HSE functional, with the experimental one. The agreement between theory and experiment is extremely good. Not only is the overall shape of the luminescence band described correctly, but all the specific features are described very accurately. In particular:

Zoom In Zoom Out Reset image size

• (i)  
  The weight of the ZPL of the theoretical spectrum ${{w}_{ZPL}}=3.8$% is very close to the experimental result ${{w}_{ZPL}}=3.2$%. Both of these quantities have been determined directly from luminescence lineshapes shown in figure 4, as discussed in section 2.2. The theoretical HR factor S = 3.67 is thus also very close to the experimental HR factor $S=3.85\pm 0.05$; the latter has been extracted from the experimental spectrum as described in section 2.2.
• (ii)  
  Both the experimental and the theoretical band show about four increasingly broad phonon replicas. The theoretical phonon frequency $\hbar {{\omega }_{0}}$ = 65 meV is in very good agreement with the experimental value $\hbar {{\omega }_{0}}$ = 64 meV.
• (iii)  
  The fine structure near the ZPL, which is representative of the coupling to acoustic phonons, agrees closely.

We conclude that calculations based on hybrid density functionals describe the vibrational properties and the luminescence lineshape of NV centres with a very high accuracy.

In figure 5 we present luminescence lineshapes calculated using all the three different theoretical approaches discussed in section 4.2. The experimental curve and the one that corresponds to the HSE functional are the same as those in figure 4. The lineshape calculated at the PBE level is qualitatively similar to the HSE one, but there are quantitative differences. In particular, the weights of the first two phonon replicas are larger, and the overall band is narrower. Figure 5 also shows that when bulk phonons are used instead, the calculated luminescence lineshape bears no resemblance to the experimental curve: it is much broader and has a very different fine structure. This result clearly shows that the consideration of the bulk phonon spectrum is not sufficient to understand the phonon sideband, challenging the discussion of [29]. Taking into account vibrational modes of the defect system is essential.

Zoom In Zoom Out Reset image size