An ab initio effective solid-state photoluminescence by frequency constraint of cluster calculation

The creation and control of coherent single photons is an important step toward building next generation quantum technology. Photons represent an optimal controllable quantum system since they can encode quantum information that interacts minimally with the environment¹ and can be easily transported and manipulated with free space and integrated photonic circuits.² Traditional single photon sources based on spontaneous parametric downconversion are probabilistic in nature and, therefore, must be heralded or post-selected, severely limiting the scalability of photonic quantum technology.³ Solid-state systems, in particular, defect centers in nanomaterials, are a promising architecture for deterministic single photon sources.⁴ These sources, however, are generally limited by decoherence due to large coupling to vibrational states, which remains one of the biggest challenges in generating and collecting high quality photon emission from defect centers in nanomaterials.⁵ The existence and frequency distribution of photon emission is observed experimentally in the photoluminescence (PL) spectrum, which is, therefore, used to investigate existing and discover new potential emitters. Ab initio calculation of PL can help to predict and explain electronic and vibrational properties of emitters. However, simulating PL for nanoscale systems is computationally hard due to the sheer number of parameters required to describe such systems, which limits the accuracy of results and the size of the systems that can be studied.

Traditionally, only electronic calculations were used for identifying and characterizing emitters.^6–8 Recently, the first ab initio calculation of the vibrational line shape was calculated for a solid-state emitter using the supercell method,^9,10 opening up the ability to identify emitters with the entire vibrationally resolved PL.¹¹ To get excited state properties, the method used $Δ$-SCF, which is only mathematically proven for excited states with different symmetry than the ground state, or with the use of an orbital-specific exchange-correlation functional.¹² While this applies to their test case of the nitrogen vacancy (NV$−$⁠) center in diamond, this does not hold for all emitters, such as the neutral nickel substitution in diamond Ni$s0$⁠.¹³ It is preferable to treat excited states with Time-Dependent Density Functional Theory (TD-DFT) as it gives formally exact excited state properties under similar approximations to DFT.¹⁴ TD-DFT has been previously used to determine the electronic spectrum of NV$−$ centers in nanodiamond^15,16 and the absorption spectrum for SiV.¹⁷ 

Here, we propose to calculate a PL spectrum for nanocrystals using the accuracy of electronic properties under TD-DFT with vibrational properties under DFT. To recover the PL spectrum for solid-state crystals, the supercell method cannot be used because it is not compatible with nanocrystals in vacuum. We introduce a low frequency cutoff that removes the vibrational modes, which are not present in the solid state. We validate the applicability of our method by calculating the PL spectrum of the NV$−$ center since it is the most studied defect in the literature. Our method can be applied to general defect systems, and it is particularly attractive for systems that cannot be solved with previous methods, such as the systems with the same symmetry in the ground and excited states.¹² 

The NV$−$ center consists of a substituted nitrogen and a removed carbon from the diamond lattice. It exhibits C$3v$ symmetry in the ground state and C$s$ symmetry in the excited state. The defect is shown in Fig. 1. The electronic and vibrational structure of the NV$−$ center has been well studied.

Electronically, the NV$−$ center consists of a triplet ground and excited state that are a linear combination of Slater determinants and cannot be simulated by DFT. Fortunately, the $ms=1$ projection can be written as a Slater determinant¹⁸ and has been shown to be sufficient for calculating transition energies.¹⁵ Configuration Interaction (CI) has been used to analyze these states.¹⁹ As CI is very computationally expensive, the TD-DFT level is a sufficient improvement on $Δ$-SCF for this work as we only investigate the triplet to triplet transition.

Vibronically, the dynamical Jahn–Teller effect of the defect has been characterized and accurate vibrational wavefunctions determined assuming a single effective frequency.²⁰ The single effective frequency assumption is implicit in the Huang–Rhys model.²¹ It has been shown²⁰ that the NV$−$ line shape can be calculated under Franck–Condon approximations by relaxing symmetry from C$3v$ to C$s$⁠. We follow that procedure for our study; therefore, it is sufficient to ignore anharmonic effects and remain in the Franck–Condon picture.

To calculate the emission spectra from this model, we assume a displaced harmonic oscillator (DHO) model under the Franck–Condon approximation. The DHO model along one normal mode is shown in Fig. 2. $Evib$ is the eigenfrequency of the mode and $D$ is the displacement of atoms on excitation along one normal co-ordinate. The vertical absorption and vertical emission (⁠$Eabs,Eemit$⁠) are the energy difference for transitions with no geometry change and are determined from TD-DFT on optimized ground and excited geometries, respectively. $EZPL$ is known as the zero-phonon-line (ZPL) in solid-state or $0→0$ transition in finite systems. $Eadiab$ is the adiabatic energy, which is the difference between the energy of the relaxed ground and relaxed excited state geometries.

In the DHO model, $EZPL$ is equal to $Eadiab$ as the ground and excited states have equal zero-point energies. $D$ for a vibrational mode defines the Partial Huang–Rhys (PHR) factor, which measures contribution of that mode to the PL. For details, see the supplementary material. The power of the displaced harmonic oscillator model is that the only parameters required to determine the PL are $Evib$⁠, $Eadiab$⁠, and the PHR factors.

Based on this, we modify the time-dependent approach to emission rate implemented by Etinski et al.²² We assume dipole emission with a rate given by Fermi’s golden rule. This rate contains an integral of a generating function over time. The Fourier transform of this function gives the PL. Our model is similar to the Huang–Rhys method used in Alkauskas et al.¹⁰ without implicitly assuming a single effective frequency.

This calculates the PL for a cluster calculation under the DHO model. The cluster calculation is affected by finite-size effects, which introduce coupling to normal modes of the surface atoms. This coupling does not exist in solid state because the crystal matrix will suppress movement along these degrees of freedom. Specifically for the solid state of $NV−$⁠, Webber et al.,²³ using the supercell method with DFT, showed that the total vibrational density of states is dramatically suppressed at low frequencies.

In order to obtain the solid-state PL spectrum from the nanocrystal calculation, we need to remove the contribution of modes at low frequencies. These can be identified from a constrained optimization that mimics solid-state conditions by fixing the position of the outermost atoms. From this, a low frequency cutoff can be identified by comparing unconstrained with constrained PHR factors. Removing the modes below the cutoff leads to a modified nanodiamond spectrum that mimics the solid-state behavior.

In summary, our method is as follows: first, a nanocrystal with defect is constructed and DFT is performed to get the ground state optimized geometries. Then, vibrational calculations under the harmonic approximation in DFT yield the $Evib$ and normal coordinates as seen in Fig. 2. The ground state is also used with TD-DFT to yield the optimized excited geometry, which gives $Eadiab$⁠, $D$⁠, and the PHR spectrum. A constrained TD-DFT geometry optimization then gives a constrained PHR spectrum, which allows us to identify the lower region of suppressed modes as well as the region of unaffected or amplified modes. The boundary between regions is the first frequencies where the constrained PHR value is less than or equal to the unconstrained. A cutoff value in this region is applied to the unconstrained PHR spectrum. The PHR spectrum with the ZPL from TD-DFT is input into the DHO model to get the PL.This method requires one ground state geometry optimization DFT calculation, one normal mode calculation, two excited state geometry optimization TD-DFT calculations with and without constraints, and the final PL calculation. These are all standard procedures in computational chemistry. The exact equations we implemented for the PL calculation are given in the supplementary material and can be implemented with minor modifications to the existing software, VIBES.²² 

Specifically, we construct a nanodiamond with composition C$197$NH$140$⁠. This corresponds to third nearest neighbors to the defect. Previous electronic structure studies^24,25 found that larger basis sets than double-$ζ$ polarization do not increase accuracy significantly for geometry optimization. For all atoms, we use def2-SV(P).²⁶ 

Gali et al.²⁴ found that PBE0 was the optimal exchange-correlation function for electronic calculations of NV$−$ under TD-DFT. As such, we use PBE0 for all calculations.

Turbomole²⁷ was used to relax the geometry using DFT²⁸ with C$3v$ symmetry constraints. All excited properties were also calculated with TD-DFT^29,30 in Turbomole using the same basis set and functional. The excited state optimization was performed under different constraints to analyze finite-size effects. First, we performed unconstrained optimization under C$s$⁠. Then, we performed optimization with all CH and CH$2$ groups constrained in ground state co-ordinates.

The program SNF³¹ was used to calculate the ground vibrational modes using a finite difference with $0.01$ Å displacements. This consists of two displacements along three degrees of freedom for $338$ atoms, which requires $2028$ data points. Redundancy due to symmetry reduces this to $427$ single shot DFT calculations.

VIBES²² was used to calculate D and subsequently the PHR factors for the transition between the ground state and both constrained and unconstrained excited states. The difference of these was used to calculate boundaries between suppressed and unaffected regions. The cutoff frequency was found as the arithmetic mean between these frequencies.

We calculate the vibrational line shape using VIBES²² modified to implement the DHO model with a cutoff. As our model includes temperature effects (see Sec. I in the supplementary material for details), the temperature was set to $300$ K. The correlation function was evaluated over $300$ fs at $16384$ intervals. To replicate experimental line broadening, the spectrum was convoluted with a Gaussian of width $200$ cm$−1$⁠, which is approximately the experimental FWHM of the ZPL in Fig. 4.

The adiabatic energy, $Eadiab$⁠, for the C$197$NH$140$ nanodiamond under TD-DFT is $1.945$ eV. This is consistent with the experimental ZPL $1.945$ eV.³³ While this appears extraordinarily accurate, it will depend on the shape of the simulated nanodiamond (see Sec. II in the supplementary material for details). Despite this, the results are an improvement on $Δ$-SCF, which reported $1.757$ eV under PBE and $2.035$ eV under HSE. These are consistent with previous TD-DFT results, which found $Eemit$ of $2.2$ eV compared to our $2.1$ eV.²⁴ 

Figure 3(a) shows the PHR factors as a function of the frequency of their respective normal mode for a C$197$NH$140$ nanodiamond. The PHR factors show the amount of coupling of a normal mode to the emission. The vibrational modes show three distinct regions. At low frequency, most modes involve large displacements from all atoms; at $480$–$650$ $cm−1$⁠, the vibrational modes have the largest amplitude with the closest atoms to the defect, as shown in Fig. 1; higher frequency modes involve smaller subregions; for instance, modes around $1500$ $cm−1$ are exclusively C–H bending modes but cannot be seen in Fig. 2 since they do not couple to emission and have a PHR factor of zero. The middle region at $532$ $cm−1$ corresponds to the $65$ meV peak identified by Alkauskas et al.¹⁰ 

In Fig. 3(b), we report the PHR factors obtained by constraining the CH and CH$2$ groups, and by plotting the difference between constrained and unconstrained factors in Fig. 3(c), we can see that at low frequencies (on the left of the blue boundary), the presence of large positive peaks indicates that the coupling to vibrational modes is present in the nanodiamond but not in bulk, whereas at high frequencies (on the right of the red boundary), the vibrational modes equally couple to the transition for both bulk and nanodiamond. The region within the boundaries, therefore, separates the unique properties of the nanodiamond from those in common with bulk. Eliminating the vibrational modes with frequencies below this cutoff region should return the bulk behavior. In Fig. S1 of the supplementary material, we compare the unconstrained and constrained PHR factors for three different nanocrystal shapes and show that the suppression of the low frequency peaks occurs under constraint across all shapes.

We can also compare the PHR factors with vibrational data obtained for bulk $NV−$ as reported in Webber et al.²³ As the authors report in Fig. 2, the total density of states begins to increase at $200$ $cm−1$⁠, illustrating a lack of vibrational states at low frequencies. Comparison with Fig. 3(a) shows that the nanodiamond has large coupling to frequencies in this same region. The coupling of low frequency modes we see in the nanodiamond can, thus, be considered unique to the nanodiamond and is not present in bulk. This motivates our use of a low frequency cutoff to eliminate these in order to imitate the solid-state behavior. Note that we derive the value of our cutoff from nanodiamond calculations with constraints, not from solid-state data. Similarly, consider the bulk vibrational states projected onto nitrogen also reported in Ref. 23. These indicate vibrational modes that involve nitrogen and are likely to be important for defect transitions. These modes begin around $350$ $cm−1$ and peak around $600$ $cm−1$⁠. The peak corresponds to a similar cluster of peaks in the PHR factors around $532$ $cm−1$ in Fig. 3(a). Furthermore, the lowest frequency nitrogen modes are similar to our lower bound for the cutoff at $357$ $cm−1$ calculated from comparing unconstrained to constrained PHR factors. The constraint of outer carbons would eliminate vibrational modes that involve the edge of the nanodiamond, but not modes of the defect, in particular, nitrogen. It is unsurprising, then, that we see correspondence between our cutoff and the $350$ $cm−1$ edge of nitrogen projected modes in bulk. The two distinct peaks of Fig. 3(a) demonstrate why the assumption of a single effective frequency from the Huang–Rhys model is invalid. In comparison, PL for solid state NV$−$ has been reproduced under the mean frequency assumption.¹⁰ In that sense, our low frequency cutoff reproduces solid-state spectra from cluster calculations by mimicking a single effective frequency around $532$ $cm−1$⁠.

The predicted PL spectra are given in Fig. 4. [Figure S2(b) in the supplementary material extends the results for a larger range of cutoff frequencies.] The spectra with no cutoff is the ab initio spectra for a C$197$NH$140$ nanodiamond in vacuum based on our displaced harmonic oscillator model. Experimental spectra for small size nanodiamonds is not available because NV$−$centers in such small crystals are unlikely to form.^34,35 In comparison, both PL spectra with cutoff match the peaks of the solid-state experimental spectrum and recreate the line shape of the sideband well for both $357$ and $469$ $cm−1$ cutoffs. The spectrum without cutoff has a larger intensity at long wavelengths than the experiment. The spectrum based on cutoff at $357$ $cm−1$ has a longer tail than the one at $469$ $cm−1$⁠, which shows that the tail is due to the contribution from the large low energy peak in Fig. 3(a) as well as from the smaller peak around $400$ $cm−1$⁠. The cutoff spectra are less intense at long wavelengths but is still higher than the experiment.

We are not able to recreate the long wavelength tail of the spectrum. This is likely due to the inelegance of using a strict cutoff on low frequency modes. It is likely that higher frequency modes must be modified. This idea is consistent with Fig. 3 where the vibrational mode around $645$ $cm−1$ increases in the constrained case. Furthermore, the constrained calculation is not an exact simulation of the solid state. We freeze the CH and CH$2$ groups when, in reality, they are only suppressed by a diamond lattice. These can be explored in future work. Nevertheless, we have established that clusters introduce a region of low frequency vibrations that couple to excitation and cause the PL to have a much higher intensity at long wavelengths compared to solid-state PL.

We have demonstrated an ab initio method for calculating the photoluminescence spectrum of solid-state emitters from cluster calculations. We have shown that finite-size effects from cluster calculations appear as part of a low frequency peak in the Huang–Rhys spectrum. We have shown that a low frequency cutoff recovers the solid-state PL spectrum. Furthermore, we have demonstrated the first vibrationally resolved PL spectrum for a nanodiamond as well as the first spectrum utilizing TD-DFT for excited state properties, which has fewer restrictions than $Δ$-SCF. This method was demonstrated with NV$−$ and found to match experimental PL peaks but loses intensity from the ZPL to long wavelengths. A more general approach to predict the value of the cutoff will require future work involving comparisons across multiple defects.

Our method enables the calculation of the photoluminescence spectrum for defects of arbitrary symmetry between the ground and excited states beyond the NV$−$ center. This constitutes a valuable tool for understanding and predicting solid-state single photon emitter properties.

See the supplementary material for an analytic derivation of the PL spectrum, electronic and vibrational results for C$121$NH$100$ and C$145$NH$100$⁠, and spectra generated by varying the cutoff frequency as well as from constrained data without cutoff.

A.P. acknowledges the RMIT University Vice-Chancellor’s Senior Research Fellowship, the ARC DECRA Fellowship (No. DE140101700), and the Google Faculty Research Award. This work was supported by the Australian Government through the Australian Research Council under the Centre of Excellence scheme (Nos. CE170100012 and CE170100026). It was also supported by computational resources provided by the Australian Government through the National Computational Infrastructure Facility and the Pawsey Supercomputer Centre.

The data that support the findings of this study are available from the corresponding author upon reasonable request.